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A  complete  continuum  thermoelastic  theory  for  large  deformation  of  crystals  of  arbitrary 
symmetry  is  developed.  The  theory  incorporates  as  a  fundamental  state  variable  in  the 
thermodynamic  potentials  what  is  termed  an  Eulerian  strain  tensor  (in  material 
coordinates)  constructed  from  the  inverse  of  the  deformation  gradient.  Thermodynamic 
identities  and  relationships  among  Eulerian  and  the  usual  Lagrangian  material  coefficients 
are  derived,  significantly  extending  previous  literature  that  focused  on  materials  with 
cubic  or  hexagonal  symmetry  and  hydrostatic  loading  conditions.  Analytical  solutions  for 
homogeneous  deformations  of  ideal  cubic  crystals  are  studied  over  a  prescribed  range  of 
elastic  coefficients;  stress  states  and  intrinsic  stability  measures  are  compared.  For 
realistic  coefficients,  Eulerian  theory  is  shown  to  predict  more  physically  realistic  behavior 
than  Lagrangian  theory  under  large  compression  and  shear.  Analytical  solutions  for  shock 
compression  of  anisotropic  single  crystals  are  derived  for  internal  energy  functions  quartic 
in  Lagrangian  or  Eulerian  strain  and  linear  in  entropy;  results  are  analyzed  for  quartz, 
sapphire,  and  diamond.  When  elastic  constants  of  up  to  order  four  are  included,  both 
Lagrangian  and  Eulerian  theories  are  capable  of  matching  Hugoniot  data.  When  only  the 
second-order  elastic  constant  is  known,  an  alternative  theory  incorporating  a  mixed 
Eulerian-Lagrangian  strain  tensor  provides  a  reasonable  approximation  of  experimental 
data. 
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1.  Introduction 

Nonlinear  continuum  thermoelasticity  provides  a  physical  description  of  behavior  of  crystalline  solids  in  the  study  of 
acoustic  and  shock  waves,  ballistic  impact,  and  high-pressure  geophysics  problems.  For  extremely  high  pressure  events, 
much  work  has  focused  on  development  of  scalar  equations  of  state,  e.g.,  scalar  relations  among  pressure,  volume,  entropy, 
and  temperature  and  associated  thermodynamic  (energy)  potentials.  For  ductile  substances  (e.g.,  many  metals)  and  those 
that  fracture  easily,  such  a  description  is  sufficient  in  many  cases  wherein  the  deviatoric  stress  (i.e.,  shear  components)  are 
restricted  in  magnitude  due  to  plastic  slip,  twinning,  or  fracture  that  limit  shear  strength  to  a  small  fraction  of  the  applied 
pressure.  However,  some  crystalline  materials  such  as  ceramics  and  hard  minerals  may  retain  significant  shear  strength  at 
finite  strain  under  high  pressure  loading,  e.g.,  as  observed  in  plate  impact  or  explosive  loading  (Wackerle,  1962;  Fowles, 
1967;  Graham  and  Brooks,  1971 ;  Kondo  and  Ahrens,  1983;  Lang  and  Gupta,  2010).  In  these  cases,  a  three-dimensional  tensor 
theory  of  nonlinear  thermoelasticity  is  required.  Furthermore,  apart  from  exceptional  cases  such  as  hydrostatic  loading  of 
cubic  crystals,  material  anisotropy  must  be  addressed  in  descriptions  of  single  crystal  behavior. 
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The  present  work  distinguishes  among  theories  based  on  what  are  labeled  Lagrangian  and  Eulerian  finite  strain 
measures.  Let  x  denote  the  time-dependent  spatial  position  of  a  material  element  located  at  point  X  in  the  undeformed 
body.  The  deformation  gradient  is  F=V0x,  where  V0  is  the  material  gradient  operator  and  x  =  x(X,t).  The  inverse 
deformation  gradient  is  F_1  =  VX,  with  V  the  spatial  gradient  and  the  inverse  motion  X  =  X(x,  t ).  The  ratio  of  current  to 
initial  volume  of  the  element  is  J  =V/V0  =  detF\  inverting  this  gives  J-1  =  V0/V  =  detF-1.  A  theory  whose  independent 
state  variable  entering  the  thermodynamic  potentials  is  constructed  from  stretch  raised  to  some  positive  power/exponent  is 
labeled  here  as  “Lagrangian”;  a  theory  utilizing  stretch  raised  to  some  negative  power/exponent  is  labeled  here  as 
“Eulerian”.  This  terminology  will  be  explained  in  more  detail  later  by  example. 

Let  a  superscript  T  denote  transposition.  Conventional  nonlinear  elasticity  for  crystalline  solids  (Wallace,  1972;  Thurston, 
1974;  Clayton,  2011a)  incorporates  the  Lagrangian  strain  measure  E(X,t)=  2(FTF-1),  often  called  the  Green-St.  Venant 
tensor  or  simply  the  Green  strain,  in  the  thermodynamic  potentials.  This  approach,  when  elastic  constants  of  up  to  third 
order  are  included,  has  been  successful  for  modeling  many  crystalline  solids  under  compression  V/V0>0.95,  including 
ceramics  and  pure  minerals  (Winey  and  Gupta,  2004;  Clayton,  2009,  2011b;  Foulk  and  Vogler,  2010;  Clayton  et  al„  2012), 
metals  (Clayton,  2005a, b,  2006;  Vogler  and  Clayton,  2008),  and  locally  heterogeneous  geologic  materials  (Clayton,  2008, 
2010a),  but  its  accuracy  degrades  at  larger  compressions  (smaller  volume  ratios)  possible  in  shock  loading  or  ballistic  events. 
In  such  cases,  elastic  constants  of  order  four  and  higher,  difficult  to  measure  and  unknown  for  most  anisotropic  crystals,  may 
be  needed  (Thurston,  1974). 

For  hydrostatic  compression  of  cubic  crystals  or  isotropic  polycrystalline  solids,  it  has  been  shown  (Birch,  1978;  Jeanloz, 
1989)  that  pressure-volume  equations  of  state  incorporating  Eulerian  volumetric  strain  measures,  i.e.,  a  series  of  term(s) 
consisting  of  V0/V  raised  to  some  positive  exponent,  are  almost  always  more  accurate  than  those  incorporating  Lagrangian 
measures  [i.e.,  dominant  term(s)  consisting  of  V /V0  raised  to  some  positive  exponent]  when  each  representation  contains 
the  same  number  of  bulk  elastic  constants.  A  canonical  example  of  an  Eulerian  description  is  the  Birch-Murnaghan  EOS 
(Birch,  1947,  1978;  Murnaghan,  1951)  which  often  demonstrates  high  accuracy  even  when  truncated  at  second  order  and  at 
third  order  is  often  almost  indistinguishable  from  the  linear  shock  velocity-particle  velocity  relation  that  applies 
exceptionally  well  for  many  shock-compressed  solids  (Jeanloz,  1989). 

The  Birch-Murnaghan  EOS  is  by  definition  restricted  to  pressure-volume  space.  A  complete  description  for  all  stress 
states  requires  a  tensor  formulation.  For  single  crystals,  as  well  as  textured  polycrystals  and  composites,  this  formulation 
must  account  for  anisotropy.  A  mathematically  and  thermodynamically  consistent  way  to  construct  such  a  description  is  to 
assign  scalar  thermodynamic  potentials  (e.g.,  free  energy  or  internal  energy)  that  are  irreducible  functions  of  requisite 
invariants  of  an  objective  finite  strain  tensor  for  the  given  material's  symmetry.  By  conjecture,  extending  arguments  for 
Eulerian  equations  of  state  to  arbitrary  stress  states  and  anisotropic  solids,  it  is  proposed  that  thermodynamic  potentials 
incorporating  an  Eulerian  strain  measure,  as  defined  above,  will  provide  analogous  advantages  in  six-dimensional  stress- 
strain  space  as  Eulerian  equations  of  state  provide  in  pressure-volume  space.1  For  example,  if  the  analogy  holds  as 
anticipated,  Eulerian  theory  with  elastic  constants  of  up  to  order  two  might  provide  comparable  accuracy  as  Lagrangian 
theory  with  constants  of  up  to  order  three.  Higher-order  elastic  constants  are  difficult  to  measure-standard  tests  include 
wave  speed  measurements  in  stressed  crystals  (Thurston,  1974;  Thurston  et  al„  1966;  Hankey  and  Schuele,  1970)  or  costly 
shock  compression  experiments  in  multiple  directions  (Graham  and  Brooks,  1971;  Graham,  1972a,  1972b)-and  have  been 
reported  for  few  low-symmetry  materials.  Third-order  constants  can  also  be  predicted  via  first  principles  calculations  (Zhao 
et  al.,  2007).  Therefore,  any  theory  that  alleviates  the  need  for  measurements  or  atomic  calculations  of  elastic  constants 
above  a  certain  order  would  be  valuable. 

The  present  paper  develops  a  theory  that  incorporates  Eulerian  finite  strain  tensor  D(x ,  t )  =  j(1-F_1F_t),  suggested  (but 
not  implemented)  for  describing  elasticity  of  anisotropic  solids  by  Murnaghan2  and  perhaps  first  implemented  in 
calculations  (of  stressed  cubic  crystals)  by  Thomsen  (1972).  Because  D  has  components  referred  to  the  reference  coordinate 
system,  it  is  invariant  under  spatial  rotations  (Davies,  1973),  and  can  be  used  in  elastic  potentials  for  anisotropic  bodies. 
Because  D  is  symmetric  and  referred  to  material  coordinates,  functional  forms  of  thermoelastic  potentials  for  anisotropic 
materials  expressed  in  terms  of  D  will  be  the  same  as  those  in  terms  of  E  (Thomsen,  1972;  Weaver,  1976).  For  example, 
elastic  constant  tensors  of  all  orders  will  have  the  same  symmetries,  though  magnitudes  of  higher-order  constants  will  differ 
between  the  two  theories.  Transformation  formulae  can  be  derived  relating  material  constants  of  the  two  theories  (Weaver, 
1976;  Perrin  and  Delannoy,  1978),  obviating  the  need  for  additional  experiments  or  quantum  calculations  if  Lagrangian 
constants  have  already  been  obtained. 

Further  clarification  of  terms  “Lagrangian”  and  “Eulerian”  is  in  order.  In  this  paper,  a  “Lagrangian”  strain  refers  to  a  tensor 
depending  on  principal  stretch  ratios  raised  to  some  positive  exponent,  while  an  “Eulerian”  strain  refers  to  a  tensor 
depending  on  principal  stretch  ratios  raised  to  some  negative  exponent.  This  is  consistent  with  terminology  adopted  in  the 
physics  and  chemistry  literature  (Thomsen,  1972;  Weaver,  1976;  Perrin  and  Delannoy,  1978;  Davies,  1974;  Nielsen,  1986). 
According  to  this  scheme  D  is  labeled  “Eulerian”  in  the  sense  that  it  is  constructed  from  the  inverse  deformation  gradient 
F_1(x,  t)  (precisely,  right  stretch  U  raised  to  the  -2  power)  and  its  field  is  implicitly  a  function  of  spatial  coordinates  x, 


1  This  conjecture  is  later  shown  to  be  true  for  ideal  cubic  solids  with  an  ambient  pressure  derivative  of  bulk  modulus  B0m4,  but  not  necessarily  true  for 
shock  compression  of  highly  anisotropic  single  crystals. 

2  Murnaghan  (1937,  p.  257)  proposed  a  strain  energy  depending  on  deformation  measure  j  =  l-2D. 
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even  though  indices  of  D  are  referred  to  a  material  coordinate  system  (associated  with  X).  In  contrast,  according  to 
conventions  often  used  in  continuum  mechanics  literature,  a  tensor  is  said  to  be  “Lagrangian"  if  it  is  referred  to  material/ 
initial  coordinates  and  “Eulerian”  if  referred  to  spatial/current  coordinates.  According  to  this  scheme  D  would  be  Lagrangian, 
while  E  would  be  Lagrangian  in  either  scheme  (right  stretch  U  raised  to  the  +2  power,  and  expressed  in  material 
coordinates).  Almansi  strain  e(x,  t )  =  j(1-F_tF_1)  entering  Murnaghan's  (1937)  theory  for  isotropic  solids  would  be  Eulerian 
according  to  either  scheme.  The  negative  of  tensor  D  has  been  referred  to  elsewhere  as  Piola  strain  (Haupt,  2000). 

Recent  work  (Clayton  and  Bliss,  submitted  for  publication)  has  demonstrated  that  under  finite  shear,  isotropic  E-based 
Lagrangian  theory  is  prone  to  intrinsic  instability  in  terms  of  attainment  of  null  eigenvalue(s)  of  an  incremental  stiffness 
matrix  (Wang  et  al.,  1993;  Morris  and  Krenn,  2000)  with  increasing  magnitude  of  third-order  elastic  constants,  regardless  of 
their  sign.  This  is  often  referred  to  as  a  “Born  instability”  (Thomsen,  1972;  Born,  1940),  though  different  elastic  stiffness 
tensors  for  perfect  crystals  have  been  suggested  as  most  appropriate  depending  on  boundary  conditions  (Hill,  1975;  Hill  and 
Milstein,  1977).  While  certain  crystals  such  as  quartz  (Gregoryanz  et  al.,  2000)  and  boron  carbide  (Clayton,  2012,  2013)  can 
demonstrate  true  physical  instabilities,  in  a  model  such  instabilities  should  result  from  material  physics  rather  than 
pathologies  associated  with  extrapolation  of  a  strain-based  theory  to  large  deformations  outside  the  domain  for  which 
elastic  properties  have  been  measured. 

Benefits  of  using  Eulerian  strain  measures  for  nonlinear  elasticity  of  isotropic  materials  were  extolled  by  Murnaghan 
(1937)  in  the  1930s  and  were  demonstrated  for  cubic  crystals  under  hydrostatic  stress  by  Birch  (1947).  Thermal  effects  were 
considered  later  in  a  D-based  Eulerian  formulation  (Davies,  1974),  and  a  mechanical  theory  for  several  non-cubic  crystals 
incorporating  D  was  initiated  and  exercised  in  the  late  1970s  (Weaver,  1976;  Perrin  and  Delannoy,  1978).  With  the  exception 
of  the  early  works  of  Murnaghan  and  Birch,  these  Eulerian  treatments  remain  obscure,  and  theoretical  developments  and 
comparisons  with  experiment  are  limited  to  hydrostatic  pressure  loading.  Nonetheless,  the  Eulerian  framework  has 
demonstrated  superior  accuracy  over  Lagrangian  theory  for  predicting  the  hydrostatic  isothermal  response  of  a  few 
anisotropic  crystals  (Weaver,  1976;  Perrin  and  Delannoy,  1978).  Despite  such  promise,  Eulerian  D-based  thermoelasticity 
theory  has  not  been  completely  developed  for  crystals  of  arbitrary  anisotropy,  and  until  now  has  been  untested  for  general 
non-hydrostatic  stress  states.  In  the  present  work,  the  theory  is  fully  developed  and  is  applied  to  several  loading  protocols, 
including  adiabatic  uniaxial  strain  conditions  (involving  simultaneous  shear  and  compression)  characteristic  of  shock-wave 
problems  (Thurston,  1974;  Perrin  and  Delannoy-Coutris,  1983). 

The  remainder  of  this  paper  is  structured  as  follows.  Derivations  are  presented  in  parallel  for  Eulerian  and  Lagrangian 
theories  in  Section  2,  including  governing  equations,  thermodynamic  identities,  intrinsic  stability  criteria,  and  material 
coefficients.  Analytical  solutions  for  homogeneous  deformations  of  a  cubic  crystal  are  examined  in  Section  3.  These 
solutions,  which  apparently  have  not  been  given  elsewhere,  apply  for  cubic  crystals  with  fully  anisotropic  second-order 
elasticity  but  symmetrized  anharmonicity  characterized  by  a  single  third-order  constant.  A  new  solution  is  derived  in 
Section  4  for  shock  compression  of  a  single  crystal  of  arbitrary  symmetry  described  by  fourth-order  Eulerian  theory. 
Predictions  of  this  solution  are  compared  with  those  of  the  Lagrangian  solution  for  three  materials  which  remain  elastic 
under  large  uniaxial  compression:  quartz,  sapphire,  and  diamond.  Conclusions  are  given  in  Section  5.  Appendices  contain 
supporting  material  on  kinematics  and  thermomechanics.  Standard  notation  of  continuum  field  theory  is  used:  vectors  and 
tensors  are  generally  written  in  bold  italic;  when  indicial  notation  is  used,  components  of  vectors  and  tensors  are  referred  to 
a  Cartesian  frame  of  reference  and  are  written  in  plain  italic,  with  summation  over  repeated  indices. 


2.  Theory 

2.2.  Nonlinear  continuum  mechanics  of  hyperelastic  solids 

Spatial  coordinates  are  related  to  reference  coordinates  by  the  motion  x  =  x(X.  t).  The  deformation  gradient  F  and  its 
determinant  are3 

F=  V0x  (Fy  =  djXj),  7  =  detF>0.  (2.1) 

Let  P(X,t)  and  MX.  f )  denote,  respectively,  first  Piola-Kirchhoff  and  Cauchy  stress: 

P  =>F-t;  Py  =JalkFj-k'  =  J  aikdkXj.  (2.2) 

Let  ( •  )  =  d(  ■  )/dt\x  denote  the  material  time  derivative  and  u  particle  velocity.  Balance  of  linear  momentum  in  the  absence  of 
body  force  and  balance  of  angular  momentum  are 


Vo  •  p=po»; 

djPiJ  =  Po*h 

(2.3) 

t 

ii 

PijFkj  =  PkjFij ■ 

(2.4) 

See  Appendix  A  for  a  thorough  discussion  of  kinematics. 
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Let  T  denote  Helmoltz  free  energy  per  unit  initial  volume,  and  let  6  and  r/  denote  absolute  temperature  and  entropy  per 
initial  volume.  Internal  energy  density  U  obeys 


u  =  r  +  % 

The  following  usual  functional  forms  are  assumed  for  homogeneous  solids: 

(2.5) 

'P  =  !F(F,  9),  U  =  U(F,rj). 

(2.6) 

Dependence  on  F  will  be  replaced  later  by  dependence  on  symmetric  strain  measures  that  respect  rotational  invariance  of 
the  thermodynamic  potentials. 

The  local  balance  of  energy,  in  the  absence  of  scalar  heat  sources,  is 

U  =  P  :  F-V o  Q;  U  =  P^jX.-d/Qj; 

with  Q  the  referential  heat  flux  vector.  The  local  entropy  production  inequality  is 

(2.7) 

n  +  V0  •  (6T1  Q)>0;  ft/  +  d/Qj-d-1  Qjdjd>0. 

Using  (2.5)  and  (2.7)  in  (2.8), 

(2.8) 

P  :  F-)/d-if'-6T1Q  ■  V0ft>0. 

Substituting  from  the  first  of  (2.6), 

(2.9) 

(P-df'/dF)  :  F-(ij  +  dP/ddie-Q  ■  V0d>0, 

from  which  the  usual  constitutive  equations  of  hyperelasticity  can  be  deduced: 

(2.10) 

P  =  dP/dF ,  //  =  —dP/dd. 

From  (2.5),  (2.6)  and  letting  6  =  8(F,?i), 

dU  dT  dT  d6  dO  dU  dT  dO  d6 

(2.11) 

OF  OF  00  OF  1  OF’  dr/  dO  dr/  ^  d ?/ 

Then,  from  the  second  of  (2.11),  it  follows  that 

(2.12) 

P  =  dU/dF ,  0  =  dU  /  dt /. 

(2.13) 

2.2.  Lagrangian  and  Eulerian  variables 


Lagrangian  Green  strain  E(X,  t)  is  defined  as 

£  =  1(FTF-1);  Eu  =  ^d,xkdJXk-Su).  (2.14) 

Eulerian  strain  D(x ,  t)  is  defined  as 

D  =  1(1-F-1F-t);  D,j  =  ^S,j-dkX,dkXj).  (2.15) 

Considered  in  parallel  are  two  thermoelastic  formulations  more  specific  than  (2.6),  one  based  on  F  and  termed  “Lagrangian”, 
the  other  based  on  D  and  termed  “Eulerian”: 


T=T[E(F),0],  U=U[E(F),t])]-,  (2.16) 

T  =  T[D(F),0],  U=U[D(F).ri)].  (2.17) 


These  thermodynamic  potentials  are  all  invariant  under  spatial  rotations  since  both  F  and  D  are  referred  to  the  material 
coordinate  system.  First  Piola-Kirchhoff  stress  in  the  Lagrangian  description  is,  from  (2.11),  (2.13),  and  (A.8), 


=  or 

kL  dE,j  dFkL  tkJdEJL 


=  Fi 


dU 

dEn 


=  FkjSjL- 


(2.18) 


Second  Piola-Kirchhoff  stress  is 


SJL  =  0T/dEJL  =  dU/dEjL  =  FJk  PkL  =JFJkFLi’  aki. 

First  Piola-Kirchhoff  stress  in  the  Eulerian  description  is,  from  (2.11),  (2.13),  and  (A.13), 

P“  "  EjSt  -Fi'(%-2D;d^  =r«S,. 

The  Eulerian  analog  of  second  Piola-Kirchhoff  stress  is 

Sjl  =  OT /dDji  =  dU  /dDji  =  FkjFiLFjMPkM  =  Fk]FkNFiLFiMSNM  =JFijFkLaik. 


(2.19) 


(2.20) 


(2.21) 
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Cauchy  stress  becomes,  from  (2.2), 


_  r-lr  c  (W  -r'F-'F-'  ^ 

aii~J  FiKFiLWi~}  Ki  Li^i 


(2.22) 


Let  c(F,F)  =  c(E,F )  =  c(D,F)  denote  specific  heat  per  unit  reference  volume  at  constant  deformation,  where  from  (2.11): 
c  =  dU/dF  =  9(dr//dF)  =  -F(d2H/dF2).  (2.23) 

The  rate  of  internal  energy  can  be  expanded  as 


U  =  (dU/dF)  :  F  +  ( dU/dri), )  =P  :F  +  F[(6,//dF)  :  F  +  (dt]  100)6],  (2.24) 

Substituting  (2.24)  and  (2.23)  into  (2.7)  leads  to 

cF  =  9(d2'¥/dFd9 )  :  E-V0  ■  Q  (2.25) 

Defining  thermal  stress  coefficients  P(E,F)  and  (HD,  9)  as 

P  =  di]/dE  =  -d2W/dEdF,  ~p  =  dr\/6D=  -d2'F  /  dDd9\  (2.26) 

d2‘F/d9dFiJ  =  -/iy  =  -p]LFiL  =  -P^F^ ;  (2.27) 

and  using  (A.20),  (2.25)  can  be  written  as 

cF  =  -0p  :  E— V0  Q,  cF  =  —9p  :  D-V0  •  Q,  (2.28) 

Second-order  tensor  Griineisen  parameters  are  defined  as 

T  =  p/c,  r  =  p/c.  (2.29) 


The  following  Maxwell-type  equalities  can  be  derived  using  procedures  in  Thurston  (1974)  and  Clayton  (2011a): 


FT  =  ( 9/c)(dt]/dE )  =  -dS/dr/  =  -d9/6E,  (2.30) 

Ft  =  (9/c)(dr//dD)  =  -dS/dij  =  -09/dD-  (2.31) 

(F/cs)d  =  ( F/cs)(0E/0F )  =  dE(dn  =  - OF/dS ,  (2.32) 

(9/cS)d  =  ( 9/cS)(dD/d9 )  =  dD/dt]  =  -dF/dS:  (2.33) 

d  =  6E/d9  =  dii/dS,  a  =  dD/69  =  dtj/dS.  (2.34) 

Analogously  to  (2.26)  and  (2.27), 

FT  =  -d2U /dEdi],  9r  =  -d2U/dDdii ;  (2.35) 

d2U/d,idFy  =  -Ffy  =  -FTJLFiL  =  -0rKLFJ(j  F^FJm ■  (2.36) 


Defining  akL  =  dFkL/dF  at  constant  P,  thermal  expansion  coefficients  are  related  implicitly  by 


2 aij  =  akL(5uFkj  +  SLJFkI ),  2 a,j  =  akLFLm(Flk  Fj^  +  FJkFlr^).  (2.37) 

Specific  heats  per  unit  volume  at  constant  deformation  (c  =  c  =  c)  and  constant  stress  (c5  =  cs  =  cs)  obey  (Thurston,  1974; 
Clayton,  2011a) 

c=F(dri/dF)\ED  =  -F(O2,F/dF2)  =  0V/dF,  (P  =  9(d, ]/dF)\^  j.  (2.38) 

Isothermal  second-order  thermodynamic  elastic  coefficients  are 


qs  — 

dE, 


dSia 
ij 


dixF 


CWL  ~  dD 


dSia 

't 


d^W 

dDiidDi 


ijouKL 


OEydEia  ’ 

Isentropic  second-order  thermodynamic  elastic  coefficients  are 
7 ^’i  dSKL 


"  dE, 


d2U 

''  dEydEK i  ’ 


Al  _  dSKL 
^  -  d% 


dtU 

’’  dDydD/a ' 


Thermal  expansion  and  thermal  stress  coefficients  are  related  by 
P  =  (dri/FE) \e  =  (dri/dS) \s  :  (OS/dE) \s  =  a  :  c\  p,j  =  aKLCaKLIJ. 


P  =  (d>l/dD)\ll  =  (dri/dS)\l>  :  ( dS/0E)\e  =  a  :  0° ,  pIJ  =  dKLt‘> 


KLIJ- 


(2.39) 

(2.40) 

(2.41) 

(2.42) 
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Specific  heats  per  unit  reference  volume  are  related,  using  the  procedure  in  Thurston  (1974),  as 

C'-C=ea\p  =  e&:ji.  (2.43) 

Isentropic  and  isothermal  coefficients  are  related,  using  Maxwell  relations,  by 

C”  =  (<5S/<5F)|„  =  (dS/dE)\e  +  (dS/dG\E)®(d9/dE)\e  =  C°  +  (0/c)p®p,  (2.44) 


C"  =  (dS/<3D)|,  =  (dS/dD) \e  +  (dS/de\D)®(d0/dD)\g  =  C9  +  (6/c)P®~p,  (2.45) 

or  in  indicial  notation, 

C IjKL  =  &1JKL  +  (.mPvPKL,  ClJKL  =  ^1JKL  +  ^fi/C )PijPkL ■  (2.46) 

Strain  energy  density,  per  unit  reference  volume,  is  defined  as  follows  for  a  homogeneous  solid: 

W(F)  =  W[E(F)]  =  W[D(F)]  =  *F(F,  0O).  (2.47) 

where  00  is  a  fixed  reference  temperature.  When  temperature/entropy  effects  are  omitted,  second-order  elastic  coefficients 
reduce  to 


C  =  d2  W /dEdE,  C  =  d2W /dDdD. 

Tangent  modulus  A (F)  is  defined  as 

A  =  dP/dF  =  d2W  /  dFdF,  AtJkL  =  dPkL/dFij  =  d2W/dFijdFkL. 
Coefficients  A  and  C  are  related  by 


Apr  =  d(FkMSML)/dFi]  =  FjNFkMCNjML  +  dikSjL  =  FiNFkMCNjML  +  SikFj^P„L. 

Similarly,  coefficients  A  and  C  are  related  by 

A i)kL  =  d(FIk  FLmFPmSip)/dFij  =  Fqi  Flk  FjmFRmFLnFPnCQRip-Fjk  PiL~FLj  Pkj-FjmFLmFkNPiN. 


Equating  (2.50)  and  (2.51), 


c  ijkl  —  Fn  1  Ffjj  Fjj  1  FNj  FI(l  F0l  F^  Fp/  C  mnop  -f F ^  FJn‘  PnL  -Fn  1 F ^  F ^  PiL  -Fnl  FK ^  FLi‘  PmJ - Fn  1  FJnj  Fu J,  PiK. 


(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 


Relations  analogous  to  (2.52)  hold  when  either  isothermal  or  isentropic  coefficients  are  used.  The  local  linear  momentum 
balance  for  a  homogeneous  elastic  solid  in  the  absence  of  body  force  becomes,  with  Atf  =  d2W(F,9)/dFdF, 

Pa*i  =  ^ijkL  4  diXp  ~Pij  4/  d.  (2.53) 

Stress  power  per  unit  reference  volume  is,  from  (A.20), 

W  =  ( dW/dF )  F  =  P  F=S  :E  =  S  :D=Ja  d  (2.54) 

Let  SW  be  a  first-order  increment  in  strain  energy  associated  with  deformation  gradient  variation  SF.  Then  analogously 
to  (2.54), 

SW  =  (dW/dF)  :  SF  =  P  :  SF  =  S  :  SE  =  S  :  SD=]a  :  Se,  (2.55) 

where 

Se  =  F~J  SE  F-1  =F  SD  FT;  Sey  =  FRj SEKLFkk  =  FiKSDKLFJL ;  (2.56) 


SE  =  ’  [Ft5F  +  (i5F)tF],  Se  =  1[(5F)F_1  +  F“T(5F)T], 


(2.57) 


The  first  equality  in  each  of  (2.56),  or  the  second  equality  of  (2.57),  can  be  used  as  a  definition  for  spatial  increment  (5e;  the 
second  equality  in  each  of  (2.56)  is  consistent  with  transformation  formulae  between  SD  and  SE  analogous  to  (A.20).  This 
definition  for  Se  is  unique  when  incremental  deformation  gradient  SF  is  prescribed,  and  is  identical  to  that  used  widely 
elsewhere  in  the  analysis  of  internal  elastic  stability  (e.g.,  Morris  and  Krenn,  2000,  their  Eq.  (9)).  Making  the  connection 
(5W<->W  dt,  it  follows  that  &?<->ddt.  Integrated  quantity  f  d  dt  is  path  dependent  (in  contrast  to  F  or  D  that  depend 

only  on  current  values  of  F  or  F-1 ),  and  is  not  used  in  analysis  in  this  paper  or  others  (Wang  et  al„  1993;  Morris  and  Krenn, 
2000;  Hill,  1975)  dealing  with  intrinsic  stability. 

Spatial  modulus  c  is  defined  as 


Cijkl  — j  FjiFjjFkKFiLCjKL  —J  Fjj  Fjj  F Rk  F u  C  ijKi.~f7nSjk—o'jiSik—<7jkSj/—r>jkSji. 


(2.58) 


The  fully  symmetric  form  of  incremental  tangent  modulus  B  is  (Wang  et  al„  1993;  Morris  and  Krenn,  2000;  Clayton,  2012) 
Bp;  =  C ijki  4-  failcSji  +  auSjh  +  cr y/<5,k  +  (7jkSu-(TijSki-okiSij) 

—J  Fp  Fjj  Fpp Fu  CijKL—^ffUcSji  -f-  t>ijSjk  -{-  (TjiSjk  -{-  ctjkSji  -{-  <TjjSki  -t-  a^Sy). 


(2.59) 
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Consider  an  elastic  solid  undergoing  homogeneous  deformation,  omitting  thermal  effects.  Intrinsic  mechanical  stability  can 
be  defined  as  local  convexity  of  strain  energy  with  respect  to  a  given  strain  or  deformation  measure  (Hill,  1975;  Parry,  1978; 
Clayton,  2013)  and  therefore  depends  on  the  choice  of  conjugate  stress-strain  variables.  For  intrinsic  stability  consistent 
with  classical  stability  under  all-around  dead  loading  (Hill,  1975), 

SP  :  SF  =  SF  :  ( d2W/dFdF )  :  SF  =  SF  :  A  :  SF  >  0odet[A]  >  0.  (2.60) 

For  intrinsic  stability  consistent  with  classical  stability  under  controlled  Cauchy  stress  (Wang  et  al„  1993;  Morris  and  Krenn, 
2000), 

Sa  :  Se  =  [( da/dE )  :  SE]  :  Se  =  Se  :  B  :  Se  >  0<=>det[B]  >  0.  (2.61) 

The  first  two  equalities  in  (2.61)  strictly  apply  only  when  F  is  symmetric,  but  the  inequalities  apply  regardless  since  since 
6x6  matrix  [B]  is  rotationally  invariant.  Symmetric  coefficients  B ,jk|  (e.g.,  evaluated  from  isentropic  thermodynamic  moduli) 
enter  the  linear  momentum  equation  for  propagation  of  small  amplitude  elastic  waves  from  a  hydrostatically  stressed  initial 
configuration  (Thomsen,  1972);  when  the  initial  configuration  is  stressed  anisotropically,  different  tangent  moduli  may  enter 
the  wave  equation  (Thurston,  1965). 


2.3.  Thermoelastic  potentials  and  material  constants 


An  unstrained  reference  state  is  defined  by  (E,0)  =  (O,0o),  (D,0)  =  (O,0o),  and  temperature  change  from  this  reference 
state  is  A 6  =  6-60.  In  what  follows,  Greek  subscripts  denote  Voigt  notation  for  symmetric  indices,  e.g.,  (■)y  =  (  )jI  <->(■)„  : 

11<->1,  22  <->2,  33  <->3,  23  =  32<->4,  13  =  31^5,  12  =  21^6.  (2.62) 

Following  the  standard  convention  (Brugger,  1964;  Thurston,  1974;  Clayton,  2011a),  shear  strain  components  contain  a  factor 
of  two,  but  stresses  and  stiffness  coefficients  do  not.  For  example,  E6  =  2EU,  DB  =  2 Du,  Se  =  S12,  C45  =  C2313.  First  consider 
free  energy  per  unit  reference  volume,  which  can  be  expressed  as  in  either  of  the  following  series  expansions  about  energy 
‘F 0  from  the  reference  state: 

Y(E,  6)  =  ‘f'o  +  C oaEa  +  ^OapEaE/i  +  .y  + - fiOaEaA0—  ^pOaijEaEfiA0—^p'oaEa(A0)2 - +  g(0),  (2.63) 


W(D,9)  =  F0  +  CS0aDa  +  ^  Co apDaDp  +  lco afirDaDfiDr  +  ■■.-pOaDaA0-±pOa//DaD/!A0-p'OaDa(A0)2—.+g(0).  (2.64) 

Letting  ( •  )lo  =  ( •  )l£  =  D  =  o,fl  =  0O’  material  coefficients  with  zero  subscripts  are  constants  evaluated  at  the  reference  state, 
which  is  assumed  stress  free: 

ro  =  W(O,0o),  CeOa  =  (dW/dEa)lo  =  0;  (2.65) 


Q 6  — 
C-Oa/l  ~ 


Pda  —  ' 


d^-T 


dEadE,i 


rf  — 

'“'fVrtfv  - 


d6xF 


°“*  l  dE„dEndEy 


d9dE,. 


Pa  all  — ' 


d9dEadEfj 


>  P'oa  =  ■ 


dPW 


d6 l2dEn 


and 


r0  =  F(0, 0O),  C0a  =  (dr/dDa)  |0  =  0; 


cs  ■ 

C-Oa/l  — 


Po  a  —  ' 


tf-t 

dDadDp 
d 2W 


d  0dD„ 


c9 

w  Oa/lr  ~ 


Poall  ~  ' 


dD,,dDj’dDj 


d9dDadD/l 


,  P’Oa  =  - 


d02dD„ 


Letting  c0  denote  a  constant  specific  heat  for  the  unstrained  material,  thermal  free  energy  is  prescribed  as 
g  =  -c09  In  (9/90)=>c0  =  -0o(d2g/d 02)\o. 


(2.66) 

(2.67) 

(2.68) 

(2.69) 

(2.70) 

(2.71) 


Internal  energy  can  be  treated  in  a  similar  way,  letting  U0  denote  internal  energy  in  the  reference  state  defined  by 
(E,  tf)  =  (0,  )/0),  (D,  if)  =  (0, i/0),  and  entropy  change  from  this  reference  state  is  Ay  =  y-y0 : 
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U(D ,  r/)  —  U o  4-  C0aDa  +  C 


Oa/l^aD/)  +  C’0a/lyD„DflD),  +  ■ 


„D(,A?/  +  —  t  OapDcDp^  +  j; /" '  oaDa(Aij)2  + - h(/]) 


(2.73) 


Material  coefficients  evaluated  at  the  unstressed  reference  state  are 
1/0  =  17(0,170),  Cga  =  (dU/3Ea)\0  =  0; 


(2.74) 


c  — 

ap  ~ 


f  cfU 

) 

r"  (  ^ 

ydEadEp 

Jo’ 

°al,r  ydEadEpdEr 

(2.75) 


f'of’oo  = 


e0r  oc/i  —  ■ 


0  dr] 

CdE„ 


o 

d3  U 


di/dEadEp 


o 

f  cPU 
0°r°a~  [dn2dE„ 


cPU 

dtjdEa 


and 


U0  =  U(0,%),  Coa  =  (d‘F/dDa)\0  =  0; 


(2.76) 


(2.77) 


c"  |  d2u 

0aP  dDadD, 


_  <PU 

0aPr  l  dDadDpdDy 


(2.78) 


dof’o  a  = 

e0r  o  „p  = 


6  drj 
CdD„ 


o 

cPU 

i dj/dDcdD/j 


36 

3D„ 


d2!} 
di/dD „ 


e0  r 


,  (  d3u 

ydti23Da 


(2.79) 


Letting  c0  denote  a  constant  specific  heat  for  the  unstrained  material,  and  noting  when  the  material  is  unstrained  that 
3h/3t]  =  0/6o, 

0oh  =  c0d0(eA"/c°-l)  =  c0(3h/dO-60)^c0  =  60[dh/(dh/dn)]\0.  (2.80) 

Expanding  the  exponential  as  a  Taylor  series  gives  the  isolated  entropic  contribution 

h=  Ai/  +  ^(Af;)2/c0  +  ^(A7/)3/cg  +  (2.81) 

Material  coefficients  defined  as  derivatives  of  either  free  or  internal  energy  with  respect  to  E  are  related  to  those  defined  as 
derivatives  of  either  free  or  internal  energy  with  respect  to  D  in  Appendix  B. 


3.  Analytical  solutions:  homogeneous  isothermal  deformation  of  a  cubic  crystal 

3.1.  Cubic  crystals 

Analytical  predictions  of  constitutive  theories  based  on  strain  measures  £  and  D  are  compared  in  what  follows.  In  Section 
3,  attention  is  restricted  to  homogeneous,  isothermal  deformation  of  solids  whose  strain  energy  functions  are  truncated  at 
third  order  in  strain,  i.e.,  (2.63)  and  (2.64)  degenerate  to 

( £-  do)  =  1,  CQ apEaEp  +  g  CQ a/]yEaE/}Ey  =  ^  C a/}EaEp  +  ^Ca^EaEpEr,  (3.1 ) 

‘EiD,  6 a )  =  -)  Cq aijDaD/i  +  g  C0 apyDaDj}DY  =  -)  C apDaDji  +  gC afiyDaDjjDY,  (3.2) 

where,  without  further  consequence,  datum  energy  f'o  =  0  has  been  assigned.  When  homogeneous  deformation  F  is 
imposed,  differences  in  predictions  of  the  two  theories  arise  due  to  differences  in  strain  measures  £(£)  and  D(F),  which  from 
(A.14)  or  (A.15)  are  second  order  in  strain,  as  well  as  differences  in  third-order  elastic  constants  indicated  in  (B.9).  Second- 

order  elastic  constants  are  equal,  as  shown  in  (B.3).  Following  the  second  equalities  in  each  of  (3.1)  and  (3.2),  6  superscripts 

and  0  subscripts  are  dropped  from  the  elastic  constants,  and  C a/i  =  C 0ap  =  C0ap- 

Cubic  crystals  have  at  most  three  independent  second-order  elastic  constants  and  can  belong  to  point  groups  falling  into 
one  of  the  two  Laue  groups;  those  belonging  to  the  Laue  group  with  higher  symmetry  have  six  independent  third-order 
constants  (Thurston,  1974;  Clayton,  2011a): 

Ci,,Ci2,  C44;  Cm,  C112,  C123,  C144,  C155,  C456. 


(3.3) 
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The  same  components  of  C„pY  and  CapY  corresponding  to  CapY  in  (3.3)  are  functionally  independent,  but  numerical  values  of 
the  same  components  of  CapY  and  CapY  can  differ.  The  latter  two  are  related  by  (B.9)  or  (B.10),  yielding 

Cm  =  Cm  +  12Cn,  Cn2  =  Cll2  +  4C12,  Ci23  =  C123, 

C144  =  C144  +  2Ci2,  C155  =  C155  +  Cn  +  C12  +  4C44,  C45S  =  C456  +  3C44.  (3.4) 

In  the  reference  state,  bulk  modulus  B0,  shear  modulus  G0,  and  Poisson  ratio  v  are  defined  as 
Bo  =  3  (Cn  +  2C12),  G0  =  2(C]i-Ci2), 

p  =  (3Bq-2Go)/(6Bo  +  2Gq)  =  jCi2/(Cn  +  Ci2).  (3.5) 


Anisotropy  ratio  A  is 

A  =  1  —  G0/C44  =  1  —  ^(Gn  — Ci2)/G44;  A  =  Ooisotropic. 


(3.6) 


Notice  that  of  the  constants  (B0,  G0,v,A),  onlyA  depends  on  C44.  Second-order  elastic  constants  can  be  expressed  in  terms  of 
(Bo.t', A)  as  follows: 

Cn  1-p  C12  _  C44  _  3  \-2v  (3  71 

Bo  \+v’  Bo  \+v'  Bo  2(1  -A)  \+v 

Requiring  the  quadratic  (in  strain)  contribution  to  energy  to  be  positive  for  all  nonzero  strains  leads  to  the  restrictions 

B0  >  0,  A<1.  (3.8) 

Combinations  of  second-  and  third-order  constants  are  related  to  pressure  derivatives  of  tangent  bulk  and  shear  moduli  at 
the  reference  state  (Thurston,  1965;  Guinan  and  Steinberg,  1974): 


B'  =dB 
0  dp 


d[’(B„  +2B12)] 


dp 


'3lbGClll  +  2Cll2+ICl23 


(3.9) 


.  dG 
0  dp 


d[i(B„-Bi2)] 


dp 


0  =  “  cd-  (Cm  —  C123  +  2Go)-l. 
6B0 


(3.10) 


Here  B ap  are  components  of  incremental  stiffness  (2.59)  in  Voigt  notation,  when  stress  is  hydrostatic  (ay  =  -pSy).  These  can 
be  converted  to  expressions  in  terms  of  C apr  using  (3.4): 


3Bn 


2  Cm  +  2C112  +  ^  C123  )  +4, 


(3.11) 


G0  =  -^-(Cm-Ci23— IIC11— Ci2)-1.  (3.12) 

t>D0 

In  some  problems  analyzed  subsequently,  certain  assumptions  are  used  to  further  reduce  the  number  of  independent 
elastic  constants.  For  a  cubic  crystal  of  the  higher  symmetry  Laue  group  also  obeying  Cauchy's  relations  (Clayton,  2011a)- 
which  in  Lagrangian  theory  correspond  to  pairwise  central  force  interactions  among  atoms  and  may  omit  thermal-kinetic 
and  zero-point  vibrational  contributions  to  stiffness-two  independent  second-order  constants  and  three  independent  third- 
order  constants  remain: 

C44  =  C)2)  C)55  =  Cll2,  C456  =  C144  =  C]23-  (3.13) 

From  (3.7),  the  first  of  (3.13)  is  equivalent  to  A  =  \-(\-2v)/(2v).  Although  (3.13)  is  incompatible  with  correspondences  (3.4), 
(3.9)  and  (3.11)  can  still  be  applied  independently.  For  an  isotropic  solid  not  necessarily  obeying  Cauchy's  relations,  two 
independent  second-order  constants  and  three  third-order  constants  also  remain: 

C44  =j(Cii-Ci2);  C144  =  2(C)i2— 0)23), 

C155  =4(Gih— C112),  C456  =  g(C1n-3C112  +  2C123).  (3.14) 

For  a  third-order  elastic  material  simultaneously  obeying  Cauchy  and  isotropic  symmetry  restrictions,  (3.13)  and  (3.14) 
applied  together,  the  response  is  parameterized  by  a  single  second-order  constant  (e.g.,  Cn  or  B0): 

C)2  =  3  CnoGo/Bo  =  5  4>v  =  4  0C11  =  |Bo;  (3.15) 

and  a  single  third-order  constant  (e.g.,  Cm  or  B0B0): 

C112  =  3C123  =  gC)n=> 

Bo  =  ~27  Gm/Bo-  G0  =  -^Cln/B0-f;  (Lagrangian) 

^0  =  —  tj  C111/B0  +  4,  G0  =  -gCm/Bo  +  |.  (Eulerian) 


(3.16) 
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Note  that  one  of  the  three  third-order  relations  in  (3.14)  becomes  redundant  when  Cauchy  conditions  (3.13)  are  applied. 
Here  and  in  what  follows4  in  Section  3,  (3.9),  (3.11),  and  (3.16)  all  apply,  but  (3.4)  no  longer  necessarily  holds  since  all 
relations  in  the  latter  are  not  generally  compatible  with  (3.16).  Behavior  analyzed  in  what  follows  in  Section  3  can  be 
described  by  the  following  energy  functions,  normalized  by  bulk  modulus  B0: 


T 

3Bq 


1  v  (4+4  +  4)  + 


(El  1^22  +  E22E33  +  E33E1I )  + 


TZaT^4  +  4  +4)-^b,  0eijEklemn 


2(1+1/)'  11  22  33'  l+i/ 

x  [8ijSKl8mn  +  Sij(SKMSLN  +  8KN8LM)  +  SKl(8im8jN  +  8IN8jM )  +  8mn(8iK8jL  +  SiL5jK)  +  8IK(8JM8LN 
+Sjl(8im8kn  +  8inSkm)  +  SiL(8jM8kn  +  8jNSkm)  +  8jK(8jM8ln  +  <5;jv<5tM)], 


■  8jn8lm ) 


(3.17) 


T 

3Bq 


1-, 


2(1  +  v) 


(Dll  +  D 22  +  D33)  +  j  +  ^(DltE>22  +  D22D33  +  D33D11)  +  1  _a  - — —  (D22  +  D 23  +  D13) 


1-/41  +1/' 


1 


—  yg (B0-4)DijDklE>mn  x  [Sij8klSmn  +  8ij(SkmSln  +  SknSlm )  +  8kl(8im8jn  +  i5/n<5/m)  +  8mn{Sik8jl  +  SuSjk) 
+Sik(Sjm8ln  +  8jn8lm)  +  8jL(Sm8KN  +  SINSKM )  +  8h(8jm8kn  +  8jn8km)  +  8jK(SjMSln  +  (5;N<5iM)]. 


(3.18) 


In  these  expressions,  normalized  quadratic  contributions  to  energy,  (l/2Bo)C0/JEc,E/j  and  (l/2B0)Ca/JDaD/J,  retain  full  cubic 
anisotropy  and  depend  on  two  dimensionless  parameters  v  and  A,  while  normalized  cubic  contributions  to  energy, 
(1/6B0)C  aprEaEpEy  and  (1/6B0)C assume  isotropic  and  Cauchy  symmetries  and  depend  only  on  dimensionless 
parameter  B0  in  (3.17)  and  (3.18). 

Values  considered  subsequently  span  a  realistic  range  for  crystalline  solids:  0<B0<8,  <!/<§,  and  -1</4<J.  When  v=\, 

the  value  B0  =  4  corresponds  to  G0  =  |,  both  of  which  are  characteristic  of  pure  polycrystalline  substances  (Guinan  and 
Steinberg,  1974;  Steinberg,  1982).  Third-order  contributions  drop  out  of  Lagrangian  and  Eulerian  energies  when  B0  =  0  and 
B0  =  4,  respectively.  The  Cauchy  relations  for  third-order  constants  are  reasonable  for  some  real  materials  such  as  noble 
metals  (copper,  silver,  gold)  wherein  closed-shell  repulsive  interactions  dominate  anharmonic  properties  (Hiki  and  Granato, 
1966).  Therefore,  the  present  model  with  property  set  (v,A,  B0)  =  (|,0,4)  would  be  a  realistic  representation  of  an 
untextured  polycrystalline  noble  metal.  Different,  stronger  Cauchy-type  relations  have  been  proposed  elsewhere  (Hiki  and 
Granato,  1966)  for  cubic  solids  wherein  the  nonlinearity  is  again  characterized  by  a  single  third-order  constant,  but  these 
other  relations  are  incompatible  with  isotropy  unless  all  third-order  constants  vanish,  so  they  could  not  be  applied  to 
describe  an  untextured  polycrystal  and  are  not  investigated  further  here. 


3.2.  Hydrostatic  loading 


For  uniform  spherical  deformation  from  initial  volume  V0  to  final  volume  V  =  JV0, 

Xi=J1/3SijXj,  dKj  =  (dJ/dFij)dKFi]=JF^dKFi]  =  0\,  (3.19) 


F=J1/31=  QQ1/31,  E=^(J2'3- 1)1,  D=  *(1-T2/3)1.  (3.20) 

Axial  components  of  strain  tensors  E  and  D  for  spherical  deformation  (3.20)  and  uniaxial  strain  (to  be  discussed  in  this 
section)  are  shown  in  Fig.  1  from  20%  expansion  to  40%  compression.  Strain  component  tends  to  become  more  strongly 
negative  relative  to  Eu  as  J  decreases. 

In  homogeneous  cubic  crystals  or  isotropic  bodies,  the  stress  state  resulting  from  such  deformation  is  hydrostatic: 

°u  =  ~Psij>  PiK  =  -pJ2/38iK,  SjK  =  -PJV38jk,  Sjk  =  -pJ5/3SJK.  (3.21 ) 

Therefore,  Cauchy  pressure  is,  using  (A.8)  and  (A.13), 


n-  1  r1/3s  -  1  r1/3  l d<F 

P  —  3J  Skk  —  3J  (  8J  dEKK 


dW  1  i-S/3  c  l,-5/3  dj 

=  ~liJ=-3J  Skk  =  ~3J  {W^Dkk 


'  F=J,/31 

Energy  densities  per  unit  reference  volume  in  (3.17)  and  (3.18)  reduce  to 

¥7/Bo=ia2/3-i)2[i4Bo(/2/3-i)]. 


F~‘  =)- 


dT 

'W 


(3.22) 


(3.23) 


■P/B0  =  I(l-r2/3)2[l-KBo-4)(l-r2/3)].  (3.24) 

These  expressions  are  independent  of  anisotropy  factor  A  and  Poisson  ratio  v\  they  also  hold  for  any  cubic  crystal,  regardless 
of  any  possible  isotropic  or  Cauchy  symmetries.  From  (3.22),  Cauchy  pressures  resulting  from  third-order  Lagrangian  and 


4  A  perhaps  more  physically  plausible,  yet  more  mathematically  cumbersome,  approach  would  impose  (3.13)  only  on  Lagrangian  constants  and  then 
use  (3.4)  to  obtain  an  alternative  set  of  Eulerian  constants. 
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Fig.  1.  Lagrangian  (En)  and  Eulerian  (Dn)  strain  components  under  spherical  and  uniaxial  deformation. 


a 


viv0 


b 


VIV0 


Fig.  2.  Normalized  (a)  energy  density  and  (b)  hydrostatic  pressure  for  spherical  deformation. 


Eulerian  models  are,  respectively, 

P/Bo  =  §  a_1/3-J1/3)[l  -|B0a2/3-l)],  (Lagrangian)  (3.25) 

p/Bo  =  |(/“7/3-J_5/3)[l  +I(B0-4)0_2/3-1)]-  (Eulerian)  (3.26) 

These  can  be  compared  with  the  well  known  respective  Birch-Murnaghan  and  Murnaghan  equations  of  state  (Birch,  1947; 
Murnaghan,  1951,  1937;  Thomsen,  1970): 

p/Bo  =  |  (j-7/3_y-5/3  )[1  +  |(b0_4)0“2/3-1)],  (Birch-Murnaghan)  (3.27) 


p/B0  =  (J  b°-1)/B0,  limp/Bo  =  -Inf.  (Murnaghan) 
b„-»o 


(3.28) 


Since  pressures  in  (3.26)  and  (3.27)  coincide,  the  Birch-Murnaghan  EOS  is  obtained  directly  from  the  present  third-order 
Eulerian  elastic  theory  based  on  strain  measure  D  when  applied  to  a  cubic  crystal  subjected  to  spherical  deformation.  Values 
of  normalized  energy  density  W/ B0  and  ‘F /B0  are  shown  in  Fig.  2(a)  for  several  values  of  B0  over  compression  range 
0.6</<l.  Normalized  pressures  from  (3.25)  to  (3.28)  are  shown  in  Fig.  2(b).  As  B0  increases,  energy  and  pressure  increase 
more  rapidly  with  decreasing  volume  for  Eulerian  theory  compared  to  Lagrangian  theory.  For  physically  characteristic  value 
B0  =  4,  Eulerian  theory  provides  much  closer  agreement  with  the  Murnaghan  EOS,  giving  a  strongly  increasing  pressure  at 
large  compression  representative  of  real  materials;  for  physically  low  value  B0  =  0,  Lagrangian  theory  nearly  coincides  with 
the  Murnaghan  EOS. 

Now  consider  incremental  modulus  B  of  (2.59),  which  for  cubic  crystals  under  spherical  deformation  reduces  to 

By kl  =  J '  ^ C ij KiSnSj] 5 kKS ii  +  p( SijSki -5ikSji -SjiSjk)  =J~7/3CijKLSiISjjSkKSiL  +  p(<%<5w  +  SikSji  +  SuSjk).  (3.29) 

This  tensor  has  the  same  symmetries  as  second-order  elastic  moduli  C  and  C,  and  thus  three  independent  components.  For 
Lagrangian  theory  with  energy  function  (3.17),  these  are 


Bn/Bo  = 


3^1/3{r^  "  ^ BoC/2/3_1)_  1  (/"2/3_1  ’  f1  ”  I Bo(/2/3_1)] }’ 


(3.30) 
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B12/B0  =  3//3|t^--Ab0(/2/3-1)  +  1c/-2/3-1) 


i-|b0c/2/3-1)]}, 


B44/Bo  =  3//3 


\—2u 


2(l-A)(l+i/)  10B°^  11  2  ^  1  11 


l-^B0C/2/3-l) 


For  Eulerian  theory  with  energy  function  (3.18),  these  are 

(B0-4)(l-r2/3)  +  ^(W2/3) 


B,2/Bo  =  3/ 
B44/B0  =  3J~7/3 


-7/3  I 

(  \-V 

9 

1 

(1  +v 

10 

-7/3  1 

r  * 

3 

1 

[i+p 

10 

-7/3  I 

r  i 

— 2v 

(B0-4)(l-r2/3)  +  ^(l-J2/3) 


1  +|(B0-4)(J-2/3-l)] 
l+|<B0-4)(T2/3-l)]}, 


2(1-A)(1  +  v) 


-A(B0-4)(l-r2/3)  +  lo-J2/3)  [l  +  |(B0-4K/-2/3-l  )j  }. 


(3.31) 

(3.32) 

(3.33) 

(3.34) 

(3.35) 


For  homogeneous  hydrostatic  loading  in  which  pressure  is  applied  incrementally,  (2.61)  is  an  exact  criterion  for  elastic 
stability  (Milstein  and  Hill,  1979),  and  there  is  no  practical  need  to  consider  other  criteria  such  as  (2.60).  Intrinsic  stability 
criterion  (2.61)  can  here  be  reduced  to  the  following  three  normalized  conditions: 


.  _  B  _  Bn  +  2B12  _  G  _  Bn-B12  _  _  B44 

B  ~  Bq  _  3B0  °’  G  “  Go  ~  2G0  °’  "  “  C44  ~ 

For  Lagrangian  theory,  left  sides  of  these  equalities  become,  explicitly, 

4b  =  (1~  B0C/2/3-l)  +  i  r2/3-l  )[1-|b0(/2/3-1)])J1/3, 

-V  -  }jv». 


(3.36) 

(3.37) 

(3.38) 

(3.39) 


Notice  AB  depends  only  on  (/,B0),  whereas  AG  depends  on  (J,i/,  B0)  and  depends  on  (j.A.v,  B0).  Also,  A,,->AG  as  A->  0. 
Analogously,  for  Eulerian  theory, 


AB  =  |i-|(B0-4)(i-r2/3)+|(i-;2/3)  i  +|(b0-4)c/-2/3-d  Jr7/3, 

i-|(Bo-4)(i-r2/3)]}r7/3, 


4g  =  {  l-i^(B0-4)(l -r2/3)  +  {^(W2/3) 


5(1-2  vy 


l-|(B0-4)(l-r2/3) 


-7/3 


(3.40) 

(3.41) 

(3.42) 


Table  1  shows  stable  domains  for  spherical/hydrostatic  strain,  defined  as  J  for  which  (3.37)-(3.47)  are  positive.  For  the 
range  of  properties  analyzed,  stability  predictions  of  either  theory  are  similar,  with  Lagrangian  theory  stable  to  slightly 
larger  compression  (at  B0  =  0)  and  Eulerian  to  slightly  larger  tension  (at  B0  =  8). 

3.3.  Uniaxial  strain 


For  uniform  uniaxial  strain  in  the  X4  direction, 
xi  =JXi,  x2=X2,  X3=X3;  d/(J  =  0: 


(3.43) 


Table  1 

Stable  domains  for  spherical  deformation  over  properly  ranges  0<B'0<8,  |!n<r<(,  and  -1<A<). 


Theory 

Ab  >  0 

Aq  >  0 

o 

A 

Lagrangian 

0.00  </<  1.13 

0.81  </<  1.06 

0.90  <J<  1.03 

oo  >  p/Bo  >  —0.12 

0.21  >p/Bo  >  -0.06 

0.11  >p/B0  >  -0.03 

Eulerian 

0.76  </<  1.23 

0.88  </<1.07 

0.93  < J  <  1 .04 

0.19  >p/B0  >-0.08 

0.12  >p/B0> -0.07 

0.07  >p/B0> -0.04 
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F=[FiJ]  =  [dXi/dXj]  = 

J  0  0" 
0  1  0 

,  F-1  =  [Fj]  =  [dX,/dx]]  = 

7 o  o" 
0  1  0 

0  0  1 

0  0  1 

tn 

II 

bn 

II 

!(/2- 1)  o  o' 
0  0  0 

o 

11 

T 

■ia-r2)  o  o' 
0  0  0 

O 

O 

O 

O 

O 

O 

(3.44) 


(3.45) 


Specifically  considered  here  is  deformation  along  a  (100)  direction,  wherein  both  referential  and  spatial  Cartesian  coordinate 
axes  are  aligned  parallel  to  cubic  axes  of  the  undeformed  crystal  lattice.  Energies  per  unit  volume  in  (3.17)  and  (3.18) 
reduce  to 


'f/b° = Sx-i)2 + S^2-1*3 = I^2-i)2 


(3.46) 


^b°=S^-2)2+§£ 


a-r2)3  =  gn-r2)2 


^-^(Bo-^i-r2) 


(3.47) 


Notice  that  energy,  and  hence  stresses  and  pressure,  depend  on  material  constants  B0,  B0,  and  v,  but  not  anisotropy  factor  A 
for  uniaxial  strain  along  (100).  Normalized  energy  densities  W  / B0  and  TJ B0  are  shown  in  Fig.  3(a)  for  several  values  of  u 
(with  fixed  B0  =  4)  over  the  compression  range  0.6</<l,  and  similarly  in  Fig.  3(b)  for  several  values  of  B0  with  fixed  v  =  |.  As 
v  decreases,  shear  modulus  G0  increases,  and  energy  increases  when  J  and  B0  are  held  fixed.  As  B0  increases,  energy 
increases  more  rapidly  with  decreasing  volume  for  Eulerian  compared  to  Lagrangian  theory  as  expected  considering  the 
rapid  increase  in  |Dn|  with  decreasing  volume  in  Fig.  1.  When  B0  =  0,  Eulerian  theory  predicts  unrealistic  negative  energy 
forj  <0.77. 

First  consider  Lagrangian  theory  based  on  strain  measure  £.  In  Voigt  notation,  nonzero  second  Piola-Kirchhoff  stress 
components  Sa  =  d¥ /dEa  are 


S^f/2-!) 


Cii  +  ^Cm(/2-l ) 


|b0C/2-1) 


\-v 

\+v 


9 

28 


b0(/2-D 


(3.48) 


S2  =s3  =  J(j2-1) 


Ci2+Jcn2C/2-l) 


=  |b0(J2-1) 


\+p  140 


9  -7 

b0(/2-D 


Cauchy  stresses  ai]=J  1FjNlPiM=J  1  FuFjM^LM  and  pressure  p  =  -  ^akk  are,  in  Voigt  notation, 


<71  =JSi  =  |  Bo  (J3-J) 


9  D  <,2  1  , 

1+«_28BoC/  _1) 
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(3.52) 


The  ratio  R  of  pressure  under  uniaxial  strain  to  that  under  spherical  deformation  at  the  same  volume  ratio  J  is  found  by 
dividing  (3.52)  by  (3.25): 


a 


v/v0 


b 


1.0  0.9  0.8  0.7  0.6 

VIV0 


Fig.  3.  Normalized  energy  density  under  uniaxial  compression  for  (a)  variable  v  and  (b)  variable  B'0. 


1996 


J.D.  Clayton  /  J.  Mech.  Phys.  Solids  61  (2013)  1983-2014 


n-j2){y 

R  = - 1 
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3  (/-1/3-y1/3) 

Now  consider  Euierian  theory  based  on  strain  measure  D.  In  Voigt  notation, 
i  r  1  ,13  ,  f 1 -v  9 


(3.53) 


St  =  ^-r2) 


c„+Jc„,(i-r2)  =|B0(i-r2) 


nonzero  stress  components  Sa  =  d'F/dD „ 

i^4(B0-4)(i-r2; 


s2=s3  =  ^(i-r2)[ci2+Jc112(i-r2)  =  §Bo(i-r2)  Y^-^Bod-r2) 

Cauchy  stresses  °B  F Li  ^Mj^LM  and  pressure  p  =  -3ctkk  are,  in  Voigt  notation, 

=r3si  =  |Bo(,-3-r5)[1^4(B0-4)(i-r2) . 
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(3.58) 


The  ratio  R  of  pressure  under  uniaxial  strain  to  that  under  spherical  deformation  at  the  same  volume  ratio  J  is  found  by 
dividing  (3.58)  by  (3.26): 


R  = 


r5-r3)- 

[,7„' 

} 

3  (r7/3-r5/3) 

1+|b0-4)(/-2/3-1) 

(3.59) 


Normalized  axial  components  of  true  (equivalently,  Cauchy  or  first  Piola-Kirchhoff)  stress  P  =  -Pn  = -a n,  positive  in 
compression,  are  compared  in  Fig.  4(a)  for  several  values  of  v  at  fixed  B0  =  4  over  0.6</<l,  and  similarly  in  Fig.  4(b)  for 
variable  B0  with  v  =  As  v  decreases,  shear  modulus  G0  increases,  and  compressive  stress  increases  when  J  and  B0  are  held 
fixed.  As  B0  increases,  P  increases  more  rapidly  with  decreasing  volume  for  Euierian  compared  to  Lagrangian  theory, 
providing  a  more  physically  realistic  representation  of  real  solid  behavior  for  V/V0<0.9  (Jeanloz,  1989).  Pressure  ratios  of 
(3.53)  and  (3.59)  are  compared  for  the  same  volume  ranges  and  property  sets  in  Fig.  5(a)  and  (b).  For  the  most  physically 
representative  case  B0  =  4,  Lagrangian  theory  predicts  R  <  1  for  J<0.9,  and  Euierian  predicts  R  >  1  for  /  <  1 .  In  analysis  of 
shock  physics  data,  material  shear  strength  is  often  estimated  as  the  difference  between  P  in  a  uniaxial  compression  test  and 
p  in  hydrostatic  compression  (Graham  and  Brooks,  1971;  Kondo  and  Ahrens,  1983).  The  present  analysis  shows  that  such  a 
procedure  would  underestimate  strength  for  Lagrangian  theory  (for  R<  1)  and  overestimate  strength  for  Euierian  theory 
(P>1). 

Linder  uniaxial  strain  deformation,  stiffness  coefficients  A  and  B  become  too  lengthy  to  write  down  individually  in  closed 
form,  but  can  easily  be  calculated  using  (2.50),  (2.51)  and  (2.59).  Intrinsic  stability  criteria  (2.60)  and  (2.61)  can  be  tested  by 
considering  the  following  inequalities  (Clayton,  2012)  that  are  necessary  conditions  for  stability  under  homogeneous  strain: 

det[A]  >  0 o+A  =  2mi„([A^])/2o  >  0;  det[B]  >  0oAB  =  2min([Ba/,])//lo  >  0.  (3.60) 


Here,  7min([.])  refers  to  the  minimum  eigenvalue  of  6  x  6  matrix  [■],  and  20  =  /Wi([Ca//])  is  the  minimum  eigenvalue  of  the 
second-order  elastic  stiffness  matrix  in  the  undeformed  material.  When  F  =  1,  Aa  =  AB  =  1.  Incremental  stiffness  B  in  (2.59) 


b 


v/v0 


Fig.  4.  Normalized  axial  true  stress  under  uniaxial  compression  for  (a)  variable  v  and  (b)  variable  B'0. 
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Fig.  5.  Ratio  R  of  Cauchy  pressure  under  uniaxial  compression  to  Cauchy  pressure  under  spherical  compression  for  (a)  variable  v  and  (b)  variable  B'0. 


Table  2 

Stable  domains  for  uniaxial  strain  over  property  ranges  0<B'0<8,  (!0  <v<  t,  and  1  <A<). 


Theory 

det[A]  >  0 

det[B]  >  0 

Lagrangian 

0.81  <)<  1.02 

0.90  <J  <1.02 

0.18  >P/B0> -0.02 

0.10  >P/B0> -0.02 

Eulerian 

0.94  <7  <1.03 

0.95  <J  <1.03 

0.06  >  P/B0  >-0.03 

0.05  >  P/B0  >-0.03 

already  has  full  Voigt  symmetry  and  can  be  written  immediately  as  a  6  x  6  matrix;  the  symmetrized  form  of  A  entering 
(3.60)  is  formed  from  converting  the  following  fourth-order  tensor  to  Voigt  notation  (Clayton,  2012): 

A 1JKL  =  y,^iJkL^ilSkK  +  hjikLSjjSkK  +  ^iJlK^ilSu  +  AjnKSjj3iL).  (3.61) 

Table  2  shows  stable  domains  for  uniaxial  strain,  defined  as  J  for  which  either  of  (3.60)  applies.  For  the  range  of  properties 
analyzed,  stability  predictions  of  either  theory  are  similar,  with  Lagrangian  theory  stable  to  slightly  larger  compression  and 
Eulerian  to  slightly  larger  tension.  First  instability  occurs  for  v=  \  and  A  =  - 1,  properties  corresponding  to  lowest  shear 
modulus  C44.  In  compression,  instability  occurs  for  largest  J  at  the  minimum  considered  value  of  ES0  =  0;  in  tension, 
instability  occurs  for  smallest  J  at  the  maximum  considered  value  of  B0  =  8. 


3.4.  Simple  shear 


For  uniform  simple  shear  of  magnitude  y  in  the  X]-X2  plane,  let 
xi  =Xi  +  7X2,  x2  =X2,  X3  =X2;  d^y  =  0; 
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(3.63) 


(3.64) 


In  this  case,  J  =  det  F=  loV=  V0,  and  nonzero  strain  components  are  E2=\y2  and  E6  =  2Eu  =  y  or  D^  =  -\y2  and 
D6  =  2Du  =  y.  Specifically  considered  here  is  shearing  along  a  (100)  direction  on  a  {010}  plane,  wherein  both  referential  and 
spatial  Cartesian  coordinate  axes  are  aligned  parallel  to  cubic  axes  of  the  undeformed  crystal  lattice.  Energy  densities  per 
unit  volume  in  (3.17)  and  (3.18)  reduce  to 
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Fig.  6.  Normalized  energy  density  under  simple  shear  for  (a)  variable  A,  (b)  variable  v  and  (c)  variable  B'0. 


Normalized  energy  densities  F /B0  and  ¥'/B0  are  shown  in  Fig.  6(a)  for  several  values  of  A  (with  fixed  B0  =  4  and  k  =  5),  in 
Fig.  6(b)  for  several  values  of  v  (with  fixed  B0  =  4  and  A=0),  and  in  Fig.  6(c)  for  several  values  of  B0  (with  fixed  v=  \  and 
A  =  0),  over  the  shear  strain  range  0<z<0.8.  As  shown  in  Fig.  6(a)  and  (b),  when  B0  =  4,  as  y  becomes  large,  energy  may 
decrease  and  even  become  unrealistically  negative  for  Lagrangian  theory.  As  shown  in  Fig.  6(c),  under  simple  shear,  F  with 
B0  =  0  is  equal  to  F  with  B0  =  4,  and  W  with  B0  =  4  is  equal  to  F  with  B0  =  0,  as  can  be  verified  by  inspection  of  (3.65)  and 
(3.66).  For  r>0.45,  F  can  become  unrealistically  negative  (indicating  unstable  behavior)  in  Lagrangian  theory  for  B0>4  and 
in  Eulerian  theory  for  B0  =  0. 

First  consider  Lagrangian  theory  based  on  strain  measure  E.  In  Voigt  notation,  nonzero  second  Piola-Kirchhoff  stress 
components  S„  =  r)F /dEa  are 

5,  .  l(C,2  +  CW  +  lc„2,4.  (3.67) 

S2  =  i(C„+C,!!),2  +  lc„„-|.5Bo,!(1l:f!;;-AB[l_^B0|.2).  (3.68) 
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Cauchy  stresses,  including  shear  stress  r  and  pressure  p,  are 
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When  71=0,  a  universal  relation  (Eringen,  1962)  for  isotropic  hyperelasticity  can  be  verified: 
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Now  consider  Eulerian  theory  based  on  strain  measure  D.  In  Voigt  notation,  nonzero  stress  components  Sa  =  d'F /dDa  are 

(3.77) 
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Cauchy  stresses,  including  shear  stress  and  pressure,  are  Cauchy  stresses,  including  shear  stress  r  and  pressure  p,  are 
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When  71=0,  a  universal  relation  (Eringen,  1962)  for  isotropic  hyperelasticity  can  be  verified: 
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For  both  Lagrangian  and  Eulerian  theories,  t  is  0(y)  and  normal  stresses  are  0 (y2). 

Normalized  shear  stress  r/B0  for  Lagrangian  and  Eulerian  theories  [(3.74)  and  (3.84)]  are  shown  in  Fig.  7(a)  for  several 
values  of  71  (with  fixed  B0  =  4  and  v=  J),  in  Fig.  7(b)  for  several  values  of  v  (with  fixed  B0  =  4  and  71=0),  and  in  Fig.  7(c)  for 
several  values  of  B0  (with  fixed  v  =  \  and  71=0),  over  the  shear  strain  range  0<y<0.8.  As  shown  in  Fig.  7(a)  and  (b),  when 
B0  =  4,  as  y  becomes  large,  shear  stress  may  decrease  and  even  become  unrealistically  negative  for  Lagrangian  theory.  As 
shown  in  Fig.  6(c),  under  simple  shear,  t  of  Lagrangian  theory  with  B0  =  0  is  equal  to  t  of  Eulerian  theory  with  B0  =  4,  and  r 
of  Lagrangian  theory  with  B0  =  4  is  equal  to  r  of  Eulerian  theory  with  B0  =  0. 

Under  shear  deformation,  stiffness  coefficients  A  and  B  can  be  calculated  using  (2.50),  (2.51 ),  and  (2.59).  Intrinsic  stability 
criteria  (2.60)  and  (2.61)  can  then  be  tested  by  considering  the  necessary  conditions 

det[A]  >  0oAa  =  2min([Aa/j])//lo  >  0;  det  [B]  >  0oAB  =  ImindB^]) //to  >  0;  (3.87) 

where  Amin([-])  is  the  minimum  eigenvalue  of  6  x  6  matrix  [•],  and  A0  =  /lmin([C„^]).  The  symmetric  form  of  A  in  (3.87)  is 
computed  with  (3.61).  Because  (2.60)  and  (2.61),  when  used  with  fully  symmetric  6x6  matrices  [A^j]  and  [Ba/J],  do  not 
strictly  apply  when  the  deformation  gradient  involves  rotation  (i.e.,  when  F^FT),  in  applying  (3.87)  the  substitution  F->  U  is 
applied,  where  U  is  the  right  stretch  tensor  from  the  polar  decomposition  of  the  simple  shear  deformation  in  (3.63)  (Clayton 
and  Bliss,  submitted  for  publication): 
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(3.88) 


With  this  substitution,  strain  energies  and  stresses  referred  to  referential  coordinates  are  unchanged,  but  Cauchy  stress 
components  may  differ  from  those  presented  analytically  above  in  the  absence  of  rotation.  Table  3  shows  stable  domains  for 
shear  strain,  defined  as  y  for  which  either  of  (3.87)  applies.  First  instability  occurs  for  i/=  5  and  A  =  -l,  properties 
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Fig.  7.  Normalized  shear  stress  for  (a)  variable  A,  (b)  variable  v  and  (c)  variable  B'0. 


Table  3 

Stable  domains  for  uniaxial  strain  over  property  ranges  0<B'0<8,  and  -1<A<§. 


Theory 

det[A]  >  0 

det[B]  >  0 

Lagrangian 

Y  <  0.05 

Y  <  0.05 

r/B0  <4.8  x  1CT3 

r/B0  <4.8  x  10~3 

Eulerian 

Y  <  0.09 

/  <  0.09 

r/B0  <7.8  x  Id"3 

r/B0  <7.8  x  10“3 

corresponding  to  lowest  shear  modulus  C44.  Eulerian  theory  is  intrinsically  stable  to  j-rsO.09  with  either  criteria,  compared  to 
the  significantly  smaller  stable  domain  of  Lagrangian  theory  (j'<0.05). 


3.5.  Summary  of  analytical  results  and  discussion 

Key  aspects  of  results  in  Sections  3.2-3.4  are  summarized  as  follows.  For  hydrostatic  compression,  for  a  physically 
characteristic  value  B0ss4,  Eulerian  theory  appears  more  realistic  than  Lagrangian  theory,  giving  a  pressure  response  closer 
to  the  Murnaghan  EOS,  while  Lagrangian  theory  fails  to  predict  a  rapidly  increasing  pressure  at  very  large  compression. 
Similar  observations  have  been  made  elsewhere  (Birch,  1947,  1978;  Jeanloz,  1989).  For  uniaxial  compression,  with  B0ss4, 
Eulerian  theory  again  offers  a  more  physically  realistic  representation;  however,  when  B0  =  0,  Eulerian  theory  can  produce 
negative  strain  energy  at  large  compression.  For  simple  shear,  Eulerian  theory  is  generally  more  stable  and  provides 
physically  reasonable  behavior  (e.g.,  monotonically  increasing  energy  with  increasing  shear  strain)  for  B  >4,  while 
Lagrangian  theory  predicts  decreasing  shear  stress  and  strain  energy  at  large  shear  for  B'>4.  The  above  statements  apply 
for  an  ideal  cubic  crystal  with  highly  symmetric  anharmonic  properties.  Deviations  may  be  expected  for  highly  anisotropic 
materials,  as  shown  in  Section  4. 

This  work  is  focused  primarily  on  comparison  of  only  two  theories  based  on  two  strain  measures  E  and  D.  Lagrangian 
E-based  theory  is  considered  because,  historically,  it  is  the  most  common  measure  used  for  anisotropic  crystal  hyperelasticity. 
Eulerian  D-based  theory  is  considered  because,  as  shown  already,  it  reduces  to  the  successful  Birch-Murnaghan  EOS  under 
hydrostatic  compression,  and  because  it  has  been  demonstrated  elsewhere  to  accurately  predict  the  response  of  anisotropic 
crystals  under  pressure  (Weaver,  1976;  Perrin  and  Delannoy,  1978).  The  present  work  investigates,  for  the  first  time,  potential 
accuracy  of  D-based  theory  for  loading  conditions  involving  deviatoric  stress  in  addition  to  pressure. 
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An  infinite  number  of  possible  strain  measures  exists  for  anisotropic  hyperelasticity,  e.g.,  constructions  consisting  of 
stretch  U  raised  to  various  exponents  are  possible.  A  likely  useful  strain  energy  potential,  not  considered  in  this  work,  would 
depend  on  the  logarithm  of  the  right  stretch,  i.e.,  W  =  W(lnU).  This  logarithmic  strain  is  of  particular  interest  since  its 
spatial  counterpart,  Hencky  strain  In  V,  is  useful  for  many  isotropic  solids  under  moderate  to  large  deformation  (Anand, 
1979).  For  isotropic  solids,  the  simple  constitutive  equality  a  =J~AdW/d  In  V  holds,  though  1/  is  not  an  appropriate  state 
variable  for  anisotropic  solids.  In  contrast,  W  =  W(ln  U )  is  appropriate  for  anisotropic  hyperelasticity,  and  though  not  often 
encountered  in  the  literature,  was  considered  in  Dluzewski  (2000).  In  that  work  Dluzewski  (2000),  transformation  formulae 
among  third-order  elastic  constants  for  Green  elasticity  and  logarithmic  elasticity  were  derived,  and  it  was  shown  that  the 
latter  tend  to  be  smaller  in  magnitude  for  several  cubic  crystals,  suggesting  greater  accuracy  of  polynomial  energy  functions 
depending  on  In  U  than  E  when  truncated  at  a  certain  order.  Whether  or  not  In  (1- based  theory  is  more  accurate  than 
D-based  theory  remains  an  open  question  that  may  be  answered  in  future  work;  superiority  of  one  measure  over  another 
likely  will  depend  on  particular  material  and  loading  regime.  Since  both  In  U  and  D  increase  rapidly  in  magnitude  with 
increasing  compression,  a  polynomial  series  for  strain  energy  in  either  strain  measure  should  account  well  for  the  rapid 
increase  in  energy,  pressure,  and  stiffness  observed  in  most  solids  at  large  compression.  Computation  of  In  U,  for  example 
via  usual  algorithms  requiring  matrix  diagonalization,  is  generally  more  cumbersome  than  D,  which  requires  only  inversion 
of  F.  Computation  of  derivatives  such  as  d  In  U /dF  needed  for  stress,  tangent  stiffness,  and  various  thermodynamic 
identities,  is  possible  but  generally  tedious  (Jog,  2009). 

4.  Shock  compression  of  low-symmetry  crystals 

Considered  in  what  follows  next  is  the  material  response  under  loading  by  an  ideal  planar  shock  wave.  Crystals  with 
homogeneous  properties  but  of  arbitrary  anisotropy  are  addressed,  i.e.,  simplifying  assumptions  made  in  Section  3.1  on 
material  symmetry  are  removed.  Generic  analytical  solutions  using  nonlinear  elastic  theories  based  on  strain  measures  E 
(Lagrangian)  and  D  (Eulerian)  are  derived  in  Section  4.1.  Specific  materials-quartz,  sapphire,  and  diamond-towards  which 
the  theories  are  applied,  and  requisite  thermoelastic  properties  are  presented  in  Section  4.2.  Results  of  the  model  as  applied 
to  these  materials  are  discussed  in  Section  4.3,  with  additional  new  developments  in  Section  4.4. 


4.1.  General  1-D  solutions 


A  shock  wave  is  represented  mathematically  as  a  propagating  surface  across  which  there  may  exist  jump  discontinuities 
in  mass  density,  particle  velocity,  strain,  stress,  entropy,  temperature,  and  internal  energy.  Considered  here  are  1-D  (i.e, 
normal  or  longitudinal)  shocks.  Quantities  associated  with  material  ahead  of  the  shock  are  labeled  with  superscript  +,  with 
material  behind  superscript  -.  Material  ahead  of  the  shock  is  assumed  to  be  at  rest,  undeformed,  unstressed,  and  at  ambient 
reference  temperature  0O.  The  jump  in  an  arbitrary  quantity  ( • )  across  the  shock  is  written  as 

[(■)]=(-r-o+-  (4.D 

In  derivations  that  follow,  the  shock  moves  at  steady  natural  velocity  in  the  X  =  X^  direction.  The  deformation  gradient  is 
uniaxial  strain  of  the  form 


F 

0 

O' 

'1  +  f 

0 

O' 

0 

1 

0 

= 

0 

1 

0 

0 

0 

1 

0 

0 

1 

F-  =  [Fij]-=  0  10=  0  1  0  ;  F+  =  l. 

Behind  the  shock,  with  x  =  x^  and  u  =  uj  the  longitudinal  particle  coordinate  and  displacement, 

„  dX  „  dU  V~  pn  ,  ,, 

F=d X=1+dX  =  1+f=J  =Vo=p=’  t=du/dX- 


(4.2) 


(4.3) 


In  the  present  work  attention  is  restricted  to  compressive  shocks,  for  which  0  <  F<1  and  -1  <  f<0,  moving  with  positive 
velocity  T>  >  0.  The  only  nonzero  components  of  Lagrangian  and  Eulerian  strain  are,  respectively, 


1 


1 


E  =  E~U=  (F^-l)  =  d  1  +  ~l  ,  D  =  D~U  =  =(1-F  )  = 


1 


1 


1  — 


1 


(1+£T 


(4.4) 


Longitudinal  force  per  unit  reference  area  (or  equivalently,  current  area  under  uniaxial  strain)  behind  the  shock  is,  positive 
in  compression, 

P=-P n  =  -KFikaikT  =  ~au  ■  (4-5) 

Though  often  referred  to  as  “Hugoniot  pressure”,  stress  P  is  generally  not  equal  to  hydrostatic  pressure  p  =  -±akk  in  a  solid 
with  shear  strength. 

Let  p  =  p~  and  v  =  07  denote  mass  density  and  particle  velocity  in  the  shocked  state.  Conservation  laws  for  mass,  linear 
momentum,  and  energy-often  referred  to  as  Rankine-Hugoniot  equations-  can  be  written,  respectively,  as  (Thurston,  1974; 


2002 


J.D.  Clayton  /  J.  Mech.  Phys.  Solids  61  (2013)  1983-2014 


McQueen  et  al.,  1970) 

Po'5>  =  p{'S>-v)<*Z  =  -v/V,  (4.6) 

P  =  p0T>u^>p0T2  =  -P/iop0o2  =  -Pi,  (4.7) 


Po=  ^(iPqu2  +  [U  ])=>[[/]  =  \p0o2.  (4.8) 

From  (4.6),  requiring  1  >J  >  0  leads  to  constraints  >  v>0.  In  energy  balance  (4.8),  the  usual  adiabatic  assumption  of  null 
heat  conduction  has  been  used,  which  is  thought  appropriate  for  elastic  materials  (Thurston,  1974)  as  well  as  relatively  weak 
shocks  in  elastic-plastic  solids  (Perrin  and  Delannoy-Coutris,  1983),  but  may  not  be  valid  for  overdriven  shocks  in  elastic- 
plastic  materials  (Wallace,  1981).  The  shock  process  is  neither  isothermal  nor  isentropic;  the  entropy  inequality  can  be 
written  as  (Germain  and  Lee,  1973) 

h/VoJ^fe>M^O.  (4.9) 

Subsequent  derivations  rely  on  internal  energy-based  constitutive  models  U(£,y)  and  U(D,y)  of  (2.72)  and  (2.73). 
Derivatives  of  these  functions  with  respect  to  strain  depend  only  on  entropy  changes  Ay  from  the  reference  state  and  hence 
are  independent  of  y0  =  y+.  Furthermore,  stress  and  temperature  depend  only  on  derivatives  of  internal  energy  with  respect 
to  strain  and  entropy,  and  hence  are  independent  of  U0.  Therefore,  to  simplify  forthcoming  derivations,  let 

Uo  =  U+  =  0,  y0=y+  =  0=>[U]  =  lT  =  U,  [y]  =  ,f  =  Ay  =  y;  (4.10) 

9+ =  (dU/d?i)+ =  0O,  0-  =  (dU/d>i)~  =0=>l$}  =0~-0+ =  0-0o  =  A0.  (4.11) 

Stress  components  thermodynamically  conjugate  to  E  or  D  are  related  to  P  via 

P  =  —FySji  =  -FS  =  -(1  +  i)S,  S  =  Sn  =  dU/dEi ,  =  dU/dE-  (4.12) 


P=-J-1F^F^S,j  =  -F~3S  =  -(1  +  i)~3S,  S  =  Sn  =  dU/dDu  =  dU/dD ; 


where  all  quantities  are  evaluated  in  the  material  behind  the  shock. 

The  following  binomial  series  (Spiegel  and  Liu,  1999)  proves  useful,  where  a  is  a  non-negative  constant: 

n-u  ,  n(n-l)  2f2  ,  n(n-l)(n-2)  3  3 


(a  +/)"  =  a"  +  nan  ’/  + 


2! 


-an~T  + 


3! 


an~3f 


n(n-l)(n-2)(n-3)  4  4  n(n-l)(n-2)(n-3)(n-4)  5  5 

+  4!  u  J  +  5,  u  J  H  ■ 

Let  a=  1.  From  (4.4),  selecting  roots  corresponding  to  i  =  0  at  E  =  D  =  0, 

i  =  -1  +  (1  +  2E)1'2  =  -1  +  (1-2D)-1/2. 


From  (4.14),  with  n=\  and /  =  2 E,  the  displacement  gradient  in  the  shocked  state  is 
£  =  F-l£2+l£3-!£4  +  g£5-, 

Similarly,  letting  n  =  -  \  and  /  =  -2D, 

i  =  D  +  §D2+%D3+fD4+fD5  +  -. 

Finally,  letting  n  =  \  and  /  =  -2D, 

(1  +  i)~3  =  (1  -2D)3/2  =  1  -3D  + 1 D2  + 1 D3  + 1 D4  +  |D5  +  .... 


(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 

(4.18) 


The  above  series  are  valid  for  -1  </<  1,  which  correspond  to  0.7<V/Vo<1.7.  Using  (4.7),  (4.10),  (4.12),  and  (4.16),  the  second 
of  (4.8)  becomes,  in  terms  of  £, 

L )=-\Pi  =  \  Si(  1  +  i)  =  1  S(E  + 1  £2-l  £3  +  I  £4-|£5  +  ■■■).  (4.19) 


Using  (4.7),  (4.10),  (4.13),  (4.17),  and  (4.18),  the  second  of  (4.8)  becomes,  alternatively  in  terms  of  D, 

U=-\Pi  =  \  S«  1  +  a-3  =  l  S(P~  2  °2-z  °3-I  D4-3D5— •).  (4.20) 

Internal  energy  functions  (2.72)  and  (2.73)-using  (2.81)  and  specialized  to  the  present  uniaxial  strain  conditions  with 
(2.74),  (2.77),  and  (4.10),  and  extended  to  fourth  order  in  strain  and  second  order  in  entropy-are 

U(£,y)  =  lc„£2  +  |cln£3  +  ±Clln£4-0o  (/TEy  +  lr,£y2  +  ^rn£2y)  +  0OJ 1  +  ^y  \  ,  (4.21) 


V(D,n)=^CuD2 


1 


1 


1 


1 


+  gCmD  +  ^CnuD  -0O  (  TiDy  +  ^r^Dij  +^fuDr] 


+  &on 


(4.22) 
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Elastic  constants  in  (4.21)  and  (4.22)  are  isentropic,  and  Voigt  notation  is  used.  The  following  simplified  notation  has  been 
used  for  material  constants  referred  to  the  reference  state,  upon  consideration  of  (B.3)  and  (B.2), 

Cn  =  Con  =  Con,  A  =r01  =foi,  r'i  =r'oi  =  f'o\-  (4.23) 

Third-  and  fourth-order  isentropic  elastic  constants  with  respect  to  strain  measures  E  and  D  are  written  as 

Cm  =  C0111,  Cim  =  C01111;  Cm  =  Cg^^,  Cmi  =  Conn-  (4.24) 

From  (B.(5),  B.10),  and  a  similar  derivation  applied  to  fourth-order  elastic  constants,  higher-order  longitudinal  constants  are 
related  by 

T’ii=/’n+4r1,  Cm=Cm  +  12Cm  Cmi  =  Cmi— 18Cm— 318Cn.  (4.25) 

Longitudinal  stress  and  temperature  in  the  shocked  state  are,  for  Lagrangian  theory, 

S  =  dU/dE  =  Cn£  + 1  Cm£2  +  g Cnn£3— 0q r](T j  +  /  +  j£T '/),  (4.26) 

0  =  dU/dri  =  0„(1  4-  v/Co-rtE-r^Eri-^TuE2);  (4.27) 

and  for  Eulerian  theory 

S  =  dU/dD  =  C„D  + 1  C,„D2  + 1  ClmD3-d0'/(£i  +  fuD  +  \r\n),  (4.28) 

e  =  dU/dii  =  <90(1  +  ri/c0-riD-r'iD?i-±fuD2).  (4.29) 

Consider  Lagrangian  theory.  Substitution  of  (4.26)  into  (4.19),  with  U  =  U,  gives 

U  =  ~ j  {0o'/(C l  +  \  E\if)}E  +  j  {Cn—  do'/Ij  (£i  +  2  r'ltj)  +  r n])£2  + 1  {Cn  +  Cm  +  &o //[(£ i  +  JC'i'/l-C n]}£3 

+2?{-6C„  +3Cm  +2C1111-0o'7[¥(A  +lrlf7)-6r„]}£4 

+A  {15Cn-6Cm  +  2Clnl  +  <V?[21(£i  +  \r\n)-\5T„]}E5  +  ....  (4.30) 

Eqs.  (4.21 )  and  (4.30)  can  be  treated  as  two  equations  in  two  unknowns  U  and  and  can,  in  principle,  yield  a  solution  for 
entropy  jump  [>/]  =  y_  =  //  in  terms  of  strain,  i.e.,  *7  =  ?/(£).  For  the  strain  energy  function  U  in  (4.21)  that  is  quadratic  in 
entropy,  such  a  solution  for  y  can  most  readily  be  obtained  using  numerical  methods  when  t;  =  (V-V0)/V0  is  prescribed. 
With  y  so  obtained,  longitudinal  stresses  S  and  P  can  then  be  acquired  immediately  using  (4.26)  and  (4.12),  noting  that 
£=1  +  f=(l  +2£)’/2. 

When  U  is  a  linear  function  of  entropy,  then  a  solution  for  > /(£)  can  be  obtained  analytically  in  closed  form  (Thurston, 


1974).  In  this  simplified  case,  most  valid  for  (c^1  -r\E)r/<2,  (4.21)  and  (4.30)  reduce  to 

U  =  j  Cii£2  +JCm£3  +  24  Cmi£4-#o(C i£  +  JCnE2)*/  +  6q r/  (4.31) 

and 

U  =  —  2  {^o£ i'?}£  +  2  (Ci  i  —dof/[2  r i  +  r n])£2  +  |{Cn  +  Cm  +  6*o'/[£ l  ~E n])£3 

+^{-6C„  +3Cm  +  2C1m-0o>/[¥A-6r11]}£4  +  ^{15C11-6Cm  +2Cnll  +0o//[21r1-15r11])£5  +  ....  (4.32) 

Writing  ;/(£)  as  a  polynomial  with  constant  coefficients  a0,a^,a2, .... 

tj  =  a0  +  Qi£  +  a2E2  +  a3E3  +  a4E?  +  a5E5  +  ■■■.  (4.33) 

Substituting  (4.33)  into  (4.31 )  and  (4.32),  equating  coefficients  of  like  powers  of  £  up  to  order  5,  and  noting  that  i/0  =  >j(0)  =  0 
from  convention  (4.10), 

a0  =  a^=a2=0,  o3  =  n^o1(3Cn  +Cm),  (4.34) 

a4  =  24^o1[”6Cn  +  3Cni  +  Cim  +  £i(3Cn  +  Cm)],  (4.35) 

as  =  ^o'llSCii— 6Cm  +  2Cim  +  7T(-9Cn  +  2Cm  +  Cnii)  +  C2(3Cn  +  Cm)].  (4.36) 

Substitution  of  entropy  jump  »;(£)=  [//],  now  known  to  fifth  order  in  strain,  into  (4.26),  (4.27),  (4.12),  and  Hugoniot 

equations  (4.6)-(4.8)  then  gives  the  longitudinal  stresses,  internal  energy  jump,  particle  velocity,  shock  velocity,  and 
temperature  completely  in  terms  of  £: 

S  =  C1iE  +  jCih£2  +  (gCnn  -Ogr^a^E3 -d0(r4a4  +  r  iia3)£4-#o(£ia5  +  rna4)E5 ,  (4.37) 

P=  -(1  +  2P)1/2S,  [l/]  =  1S[(1  +  2£) — (1  +  2£)1/2],  (4.38) 

»  =  {(S/p  0)[(1  +  2£)-(l  +  2£)1(2])1/2,  (4.39) 
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©  =  {(S/A)>[(1  +  2£)-(l  +  2E)1/2]}1/2[1-(1  +  2E)1/2]-1, 


(4.40) 


6)=40(l-r1£-lr„£2).  (4.41) 

Now  consider  Eulerian  theory.  Substitution  of  (4.28)  into  (4.20),  with  U  =  U,  gives 

(>  =  4(0 o^i  +Jr'i’/))D  +  j{Cn  +  0o</[!(/’,  +ir1^)-fii]}D2  +  i{-3c11  +  Cm  +  0o>/[(/T  +ir^)  +  3f„])D3 

+1  j-6C„-9C,„  +  2C„1,  +  0O>/[|  (r,  + 1  rlf?)  +  6f  „]}D4  +  i  (-9C,,  -6Cnl-6C,m  +  0Of/[9(/T  + 1 T'l,,)  +  9f„])D5  +  .... 

(4.42) 

Eqs.  (4.22)  and  (4.42)  are  two  equations  in  two  unknowns  U  and  and  can  yield  a  solution  for  entropy  jump  ;/  =  >/(£>).  For  0  in 
(4.22)  that  is  quadratic  in  entropy,  such  a  solution  can  most  readily  be  obtained  using  numerical  methods  when  f  is  prescribed. 
Longitudinal  stresses  S  and  P  can  then  be  acquired  immediately  using  (4.22)  and  (4.42),  noting  that  F  =  1  +  £  =  (1-2D)-1/2. 


When  0  is  a  linear  function  of  entropy,  then  a  solution  for  q(D)  can  be  obtained  analytically  in  closed  form.  In  this 
simplified  case,  most  valid  for  (Cg1-T'1£)W<«2,  (4.22)  and  (4.42)  reduce  to 

0  =  2  Ci(D2  +  g  CniD3  +24  CiniD4-do(P +  JO nD2)^  +  do  1  (4.43) 

and 

0  =  4  {d0)/Pi}D  +  1  {Cii  +  d0//[|r,— fn])D2  +  |{-3Cn  +  Cm  +  d0^[r !  +  3fn])D3 

+^{-6C11-9C111+2Clm+d0»/[|A+6f11])D4  +  i{-9Cn-6C111-6C11„+dor/[9ri  +  9f  „]}D5  + -.  (4.44) 

Writing  >;(D)  as  a  polynomial  with  constant  coefficients  b0,b^,b2, ..., 

//  =  b0  +  biD  +  b2D2  +  b3D3  +  b4D4  +  b5D5  +  ■■■.  (4.45) 

Substituting  (4.45)  into  (4.43)  and  (4.44),  equating  coefficients  of  like  powers  of  D  up  to  order  5,  and  noting  that  >;0  =  r/(0)  =  0 
from  convention  (4.10), 

b0  =  bj=b2  =  0,  f>3  =  +  Cm),  (4.46) 

£>4  =  1  [~6Cn  — 9(i  in  +  Cnn  +  r  i(-9Cn  +  Cm)],  (4.47) 


4  =  ^0o1[~9Cii-6Ciii-6Ciih  +  Pi(-33Cn-6Cin  +  Cnn)  +  P2(-9Cn  +  Cm)]. 


(4.48) 


Notice  that  a3  =  b3.  Substitution  of  entropy  jump  i/(D)  =  [»/],  now  known  to  fifth  order  in  strain,  into  (4.28),  (4.29),  (4.13), 
and  Hugoniot  equations  (4.6)-(4.8)  then  gives  the  longitudinal  stresses,  internal  energy  jump,  particle  velocity,  shock 
velocity,  and  temperature  completely  in  terms  of  D: 

S  =  CnD  +  2  CinD2  +  (^C^u3-do+^b3)D3 -do(F\b4  +  t ^3b3)D4-Oo(r3b3  +  t ub4)D5,  (4.49) 


P  =  -(1-2D)3/zS,  [U]  =1S[(1-2D)-(1-2D)3/2],  (4.50) 

v  =  )(S/p0)[(l  -2D)-(1-2D)3/2])1/2,  (4.51 ) 

{(S/p0)[(1-2D)-(1-2D)3/2])1/2[1-(1-2D)-1/2]-1,  (4.52) 

d=do(l-PiD-JPnD2).  (4.53) 


From  (4.37)  and  (4.49),  contributions  to  stresses  S  and  S  from  entropy  production  are  0(£3)  and  0(D3),  respectively.  The 
foregoing  analytical  solution  for  the  elastic  shock  response  in  Lagrangian  theory  was  derived  by  Thurston  (1974);  the 
analogous  full  and  detailed  derivation  for  Eulerian  theory  has  not  appeared  elsewhere,  to  the  author's  knowledge. 

In  order  to  apply  the  above  solutions  to  particular  materials,  the  following  six  independent  constants  are  needed  at  the 
unstressed  ambient  state  at  temperature  0O:  isentropic  elastic  constants  Cm Cm, Cnn;  Gruneisen  parameters  r i ,  rn ;  and 
mass  density  p0.  Higher-order  constants  in  the  Eulerian  theory  (Cm.CnnTn)  can  be  obtained  from  those  of  Lagrangian 
theory  via  (4.25),  or  vice  versa,  or  each  set  can  be  fit  independently  to  material  data. 

In  the  application  that  follows  in  Section  4.2  in  which  Lagrangian  and  Eulerian  elasticity  theories  are  compared,  greatest 
emphasis  is  placed  on  evaluation  of  the  mechanical,  rather  than  thermal,  response,  consistent  with  internal  energy 
functions  (4.31)  and  (4.43)  quartic  in  strain  but  linear  in  entropy.  Note  from  (2.37)  that  thermal  expansion  coefficients  in  the 
reference  state  are  equivalent  in  Lagrangian  and  Eulerian  theories: 
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Perhaps  most  often  measured  is  specific  heat  at  constant  stress/pressure  cj)  =  (^(0, (90),  which  can  be  used  in  (2.43)  to  obtain 
c0  =  c(0,  d0): 

Co  —  ^0~^Oaaafl^Qa/p  ^Oa/i  =  C(^J— {@o/C())PoaPop’  (4.55) 

where  second-order  isothermal  elastic  constants  are  obtained  from  isentropic  constants  C„/;  via  (2.46).  Griineisen 
parameters  can  be  calculated  from  second-order  elastic  constants,  specific  heat,  and  thermal  expansion  coefficients 
(Thurston,  1974;  Clayton,  2011a): 

G»  =  PoJco  =  C°0apap/c0  =  Capap/c^.  (4.56) 

Experimental  measurements  of  higher-order  Griineisen  parameters  are  scarce.  A  typical  assumption  (Wallace,  1980)  for 
weak  shocks  in  crystals  is  pTijfnp0r0lj,  which  yields 

r'r,j  =  constant ^dT,j/dEKL  =  -Jr,j(dJ-'  /  dEKL)  =  F^F^Tg 

=>A  ijki  =  2  (drij/dEKL  +  dr  Ki/d£/j)lo»2(C  oij^kl  +  G oklSij).  (4.57) 

For  constant  specific  heat,  positive  Tap  correlates  decreasing  second-order  elastic  coefficients  with  increasing  temperature. 
For  a  cubic  crystal  with  scalar  Griineisen  parameter  T=lTKK,  this  assumption  corresponds  to  (d  In  T/d  In  V)lo  =  1. 
Isentropic  third-order  constants  Cap,  can  be  computed  at  the  reference  state  from  mixed  third-order  constants  C^, 
measured  in  ultrasonic  experiments  using  the  relation  (Brugger,  1964) 

Co*  =  c li  +  floGr[C^-(aC^/afl)|0].  (4.58) 

The  difference  between  isentropic  and  mixed  coefficients  is  often  smaller  than  uncertainty  in  experimental  measurements 
of  either.  In  summary,  the  five  parameters  Cn,  Cm, Cnii.AGii  can  be  calculated  for  a  material  of  arbitrary  symmetry 
using  the  above  relations  given  C^p  and  ap  (where  p=  1,2,  ...6),  q),  dC''u/dd,  C’{ , , ,  and  Cun.  Experimental  values  exist  for  a 
number  of  crystals  for  all  parameters  except  Cun;  reported  measurements  of  the  latter  are  scarce. 

4.2.  Materials 

The  theory  and  analytical  solutions  derived  in  Section  4.1  are  applied  to  analyze  shock  compression  behavior  of  single 
crystals  of  three  hard  minerals:  quartz  (ct-Si02),  sapphire  (a-Al203  or  corundum),  and  diamond  (C).  These  materials  are 
considered  because  their  ratios  of  Hugoniot  Elastic  Limit  (HEL)  to  longitudinal  elastic  moduli  are  relatively  large,  meaning 
that  elastic  deformations  in  excess  of  several  percent  volumetric  compression  can  be  achieved  in  uniaxial  compression  prior 
to  activation  of  any  inelastic  deformation  mechanisms  that  could  render  the  analysis  of  Section  4.1  physically  unrealistic.  In 
contrast,  the  nonlinear  elastic  analysis  of  Section  4.1  could  be  applied  to  more  ductile  materials  with  a  lower  HEL-e.g., 
metals  that  undergo  plastic  slip  or  deformation  twinning-but  would  be  physically  realistic  only  at  smaller  compressions 
where  effects  of  higher-order  moduli  may  be  less  evident.  Above  the  HEL,  closed-form  analytical  solutions  for  anisotropic 
solids  become  intractable  because  neither  elastic  nor  plastic  deformation  are  one-dimensional,  and  entropy  production  from 
inelasticity  can  be  substantial.  Quartz,  sapphire,  and  diamond  also  belong  to  the  limited  set  of  anisotropic  crystals  whose 
third-  and  fourth-order  elastic  constants  have  been  reported. 

Specifically,  analytical  solutions  are  compared  for  anisotropic  nonlinear  elastic  uniaxial  shock  compression  involving 
internal  energy  functions  (4.31)  and  (4.43)  (Lagrangian  and  Eulerian  theories,  respectively)  incorporating  elastic  constants 
up  to  fourth  order.  Quartz  and  sapphire  have  trigonal  (i.e„  rhombohedral)  symmetry.  Quartz  is  analyzed  for  compression 
along  the  a-axis  (X-cut,  [1210]),  b-axis  (Y-cut,  [1010])  and  c-axis  (Z-cut,  [0001]);  sapphire  is  analyzed  for  compression  along 
the  a-axis  (X-cut)  and  c-axis  (Z-cut).  Diamond  is  cubic  and  is  analyzed  for  compression  along  a  cube  axis  (X-cut,  [100]). 
Elastic  constants  are  interchanged  as  needed  for  consistency  with  notation  of  Section  4.1.  For  example,  for  c-axis  (i.e.,  Z-cut) 
uniaxial  shock  compression,  the  analysis  of  Section  4.1  remains  valid  with  Cn  replaced  by  C33,  Cm  by  C333,  r,  by  r3,  etc. 

Requisite  material  properties  are  listed  in  Table  4  corresponding  to  an  ambient  temperature  of  295  K.  Isentropic  second- 
order  elastic  constants  for  all  three  materials  are  obtained  from  experiment  (McSkimin  et  al„  1965,  1972;  Hankey  and 
Schuele,  1970),  mixed  third-order  constants  are  obtained  for  quartz  and  sapphire  from  experiment  (Thurston  et  al.,  1966; 
Hankey  and  Schuele,  1970)  and  then  converted  to  isentropic  constants  using  (4.58).  Fourth-order  Lagrangian  constants 
shown  for  quartz  and  sapphire  are  reported  from  fits  to  shock  compression  experiments  (Fowles,  1967;  Graham,  1972a)  and 
are  inherently  adiabatic.  Quartz  is  piezoelectric;  constants  listed  correspond  to  open-circuit  conditions  (i.e.,  constant  electric 
displacement).  For  diamond,  third-  and  fourth-order  constants  are  obtained  verbatim  from  quantum  mechanical 
calculations  (Nielsen,  1986);  no  attempt  is  made  to  adjust  these  for  finite  temperature.  Griineisen  parameters  are  calculated 
via  (4.56)  and  (4.57)  using  experimentally  determined  specific  heats  at  constant  pressure  (McSkimin  et  al.,  1965;  Furukawa 
et  al.,  1956;  DeSorbo,  1953),  linear  thermal  expansion  coefficients  (McSkimin  et  al.,  1965;  Burghartz  and  Schulz,  1994;  Slack 
and  Bartram,  1975),  and  isentropic  second-order  elastic  constants.  For  quartz  and  sapphire,  third-order  Eulerian  coefficients 
Cm  are  obtained  using  conversion  (4.25).  As  discussed  later  in  Section  4.3,  fourth-order  Eulerian  constants  Cnn  for  quartz 
and  sapphire  are  fit  to  shock  velocity  versus  particle  velocity  data  (Fowles,  1967;  Graham  and  Brooks,  1971)  keeping  third- 
order  elastic  constants  fixed,  following  the  same  procedure  used  for  Cun. Fitting  this  constant  independently  rather  than 
using  the  last  of  (4.25)  provides  for  the  most  fair  comparison  of  fourth-order  Lagrangian  and  Eulerian  theories.  For  diamond, 
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Table  4 

Physical  properties  of  single  crystals  (60  =  295  K;  p0  in  g/cm3;  C,rjl  in  GPa). 


Property 

Quartz  (X) 

Quartz  (Y) 

Quartz  (Z) 

Sapphire  (X) 

Sapphire  (Z) 

Diamond  (X) 

Cn 

87.6 

87.6 

106 

497 

498 

1079 

Cm 

-211 

-333 

-814 

-3870 

-3340 

-6300 

Cm 

840 

718 

455 

2090 

2640 

5570 

Cmi 

15  930 

15  930 

18  490 

50  000 

50  000 

43  600 

15  000 

10  500 

6500 

10  000 

20  000 

16  300 

n  =/’„ 

0.74 

0.74 

0.58 

1.29 

1.29 

0.81 

ru 

3.70 

3.70 

2.90 

6.43 

6.46 

4.04 

Phel/Ch 

0.10 

0.10 

0.15 

0.05 

0.05 

0.08 

Po 

2.65 

3.98 

3.51 

B’o 

6.3 

4.2 

4.0 

u/C0  o/C„ 


Fig.  8.  Predicted  and  experimental  (Fowles,  1967)  shock  velocity  versus  particle  velocity  for  quartz,  normalized  by  linear  elastic  wave  speed  C0:  (a)  X-  and 
Y-cut  and  (b)  Z-cut. 


Eulerian  higher-order  constants  Cm  and  Cmi  are  again  taken  verbatim  from  Nielsen  (1986),  where  they  have  been 
obtained  by  fitting  numerical  data  directly  rather  than  using  (4.25).  For  each  material,  pressure  derivatives  of  bulk  modulus 
B0  (Thurston  et  al„  1966;  Nielsen,  1986;  Clayton,  2009)  are  also  listed.  Finally,  maximum  HEL  stresses  PHel  from  shock 
experiments  (Wackerle,  1962;  Fowles,  1967;  Graham,  1972b;  Lang  and  Gupta,  2010)  are  shown  for  reference,  normalized  by 
second-order  moduli  for  shocks  in  corresponding  directions.  The  domain  of  validity  of  elastic  analysis  can  be  estimated  as 
V/Vo>(C,i-PHel)/Cii. 


4.3.  Results  and  discussion 

Predicted  shock  velocity  35  versus  particle  velocity  u  is  compared  with  experimental  shock  compression  data  of  Fowles 
(1967)  in  Fig.  8  for  X-,  Y-,  and  Z-cut  quartz  specimens.  Experimental  data  are  obtained  from  plane-wave  explosive  loading 
tests  in  which  two-wave  structures  were  often  recorded  (Fowles,  1967).  Data  considered  here  correspond  only  to  the  first, 
elastic  shock  wave  in  such  tests.  Velocities  are  normalized  by  longitudinal  linear  elastic  wave  speed 

Co=(C„/Vo)1/2.  (4.59) 

Lagrangian  fourth-order  constant  Cmi  for  each  orientation  (Table  4)  was  fit  to  the  data  in  Fowles  (1967).  Eulerian  fourth- 
order  constant  Cmi  in  Table  4  has  been  fit  to  this  same  data  in  an  analogous  fashion  here.  Both  Lagrangian  and  Eulerian  fits 
are  considered  adequate  for  each  orientation.  The  unusual  nonlinearity  (i.e.,  curvature)  in  the  35-d  data  was  noted  in  Fowles 
(1967);  Eulerian  theory  predicts  relationships  with  greater  curvature.  Hugoniot  stress  (i.e.,  P)  normalized  by  Cn  is  shown  for 
each  orientation  in  Fig.  9,  along  with  experimental  data  (Fowles,  1967).  Predictions  marked  “4th  order"  are  obtained  using 
complete  solutions  and  all  material  constants.  Predictions  marked  “3rd  order”  assume  Cmi  =0  or  Cmi  =  0.  Predictions 
marked  “2nd  order”  assume  Cm  =  Cmi  =0  or  Cm  =Cmi  =0-  These  designations  apply  for  respective  Lagrangian  or 
Eulerian  solutions.  Predictions  marked  “2nd  order  mixed"  are  discussed  later  in  Section  4.4.  For  each  orientation,  4th  order 
theories  are  required  to  accurately  match  the  experimental  Hugoniot  data;  2nd  and  3rd  order  models  are  insufficient. 

Predicted  shock  velocity  35  versus  particle  velocity  o  is  compared  with  experimental  shock  compression  data  of  Graham 
and  Brooks  (1971)  in  Fig.  10  for  X-  and  Z-cut  sapphire.  Experimental  data  are  obtained  from  flyer-plate  and  plane-wave 
explosive  loading  configurations;  in  the  latter,  two-wave  structures  were  sometimes  generated  (Graham  and  Brooks,  1971). 
Data  considered  here  correspond  only  to  the  elastic  shock,  with  the  secondary,  slower  “plastic”  wave  in  which  the  HEL  was 
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Fig.  9.  Predicted  and  experimental  (Fowles,  1967)  Hugoniot  stress  for  quartz,  normalized  by  longitudinal  second-order  elastic  constant  Cu:  (a)  X-cut, 
(b)  Y-cut  and  (c)  Z-cut. 


o/C0 


Fig.  10.  Predicted  and  experimental  (Graham  and  Brooks,  1971 )  shock  velocity  versus  particle  velocity  for  sapphire,  normalized  by  linear  elastic  wave  speed  C0. 
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Fig.  11.  Predicted  and  experimental  (Graham  and  Brooks,  1971 )  Hugoniot  stress  for  sapphire,  normalized  by  longitudinal  second-order  elastic  constant  Cu : 
(a)  X-cut  and  (b)  Z-cut. 
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exceeded  not  addressed.  Velocities  are  normalized  by  wave  speed  (4.59).  Lagrangian  fourth-order  constant  Cml  for  each 
orientation  (Table  4)  was  fit  to  this  experimental  data  in  Graham  (1972a).  Eulerian  fourth-order  constant  Cnn  in  Table  4  has 
been  fit  to  this  same  data  here.  Considering  scatter  in  the  data,  both  Lagrangian  and  Eulerian  fits  are  adequate  for  each 
orientation,  giving  nearly  linear  T>-v  curves  over  the  relatively  small  compression  range  for  which  sapphire  remains  elastic 
(VVVo>0.95).  Hugoniot  stress  (i.e.,  P)  normalized  by  Cn  is  shown  for  each  orientation  in  Fig.  11,  along  with  experimental 
data  (Graham  and  Brooks,  1971).  For  each  orientation,  3rd  and  4th  order  Lagrangian  and  Eulerian  theories  are  all  capable  of 
accurately  matching  the  experimental  stress  data.  Sufficiency  of  3rd  order  Lagrangian  theory  was  also  noted  in  previous 
work  (Clayton,  2009).  Second-order  elastic  models  are  inaccurate,  with  2nd  order  Eulerian  theory  too  stiff  and  2nd  order 
Lagrangian  theory  too  compliant. 

For  predictions  of  the  shock  response  of  diamond,  all  higher-order  elastic  constants  in  Table  4  have  been  taken 
directly  from  the  quantum  mechanical  results  of  Nielsen  (1986)  since  experimental  measurements  of  third-order 
constants  of  diamond  apparently  have  not  been  reported.  Predictions  of  normalized  shock  velocity  and  Hugoniot 
stress  are  given  in  Fig.  12(a)  and  (b),  respectively,  compared  with  experimental  data  of  Lang  and  Gupta  (2010).  These 
data,  obtained  from  flyer-plate  experiments,  consist  of  five  tests  for  which  a  high  HEL  was  observed  (peak  shock 
pressures  of  «90  GPa)  and  six  corresponding  to  peak  shock  pressures  in  excess  of  f»115  GPa,  which  demonstrated  a 
marked  reduction  in  HEL  strength  (Lang  and  Gupta,  2010).  No  fitting  or  adjustment  of  third-  or  fourth-order 
Lagrangian  or  Eulerian  constants  has  been  undertaken,  so  the  comparison  of  results  can  be  deemed  as  much  a 
confirmation  of  accuracy  of  atomic  calculations  as  a  test  of  relative  merits  of  Lagrangian  and  Eulerian  theories  of 
various  orders.  From  Fig.  12(a),  4th  order  Lagrangian  theory  provides  a  better  fit  to  shock  velocity  than  4th  order 
Eulerian  theory  at  larger  particle  velocities  corresponding  to  the  higher  HEL,  with  4th  order  Eulerian  theory  tending 
to  overestimate  35.  Conversely,  at  smaller  particle  velocities  corresponding  to  the  lower  HEL,  4th  order  Eulerian 


a  b 


Fig.  12.  Predicted  and  experimental  (Lang  and  Gupta,  2010):  (a)  shock  velocity  versus  particle  velocity,  normalized  by  linear  elastic  wave  speed  C0  and 
(b)  Hugoniot  stress,  normalized  by  longitudinal  second-order  elastic  constant  Cn. 


Table  5 

Thermodynamic  predictions. 


Material 

Shock  direction 

V  JVq 

e/Go 

Lagrangian 

o/e o 

Eulerian 

n/4 

Lagrangian 

Eulerian 

P/P" 

Lagrangian 

P/P" 

Eulerian 

Quartz 

X 

0.96 

1.028 

1.028 

0.002 

0.002 

1.0002 

1.0002 

0.92 

1.055 

1.052 

0.030 

0.050 

1.0013 

1.0024 

0.88 

1.079 

1.069 

0.133 

0.418 

1.0032 

1.0039 

Y 

0.96 

1.028 

1.028 

0.003 

0.002 

1.0002 

1.0002 

0.92 

1.055 

1.052 

0.036 

0.039 

1.0016 

1.0019 

0.88 

1.079 

1.069 

0.153 

0.300 

1.0036 

1.0079 

Z 

0.96 

1.022 

1.022 

0.007 

0.006 

1.0004 

1.0004 

0.92 

1.043 

1.041 

0.066 

0.071 

1.0017 

1.0017 

0.88 

1.062 

1.054 

0.246 

0.373 

1.0033 

1.0042 

Sapphire 

X 

0.96 

1.049 

1.049 

0.016 

0.015 

1.0008 

1.0007 

0.92 

1.095 

1.091 

0.141 

0.135 

1.0028 

1.0024 

0.88 

1.137 

1.119 

0.483 

0.531 

1.0053 

1.0055 

Z 

0.96 

1.050 

1.049 

0.014 

0.012 

1.0007 

1.0006 

0.92 

1.095 

1.091 

0.122 

0.121 

1.0025 

1.0024 

0.88 

1.137 

1.120 

0.430 

0.578 

1.0049 

1.0071 

Diamond 

X 

0.96 

1.031 

1.031 

0.033 

0.042 

1.0003 

1.0003 

0.92 

1.060 

1.057 

0.255 

0.311 

1.0009 

1.0010 

0.88 

1.086 

1.075 

0.811 

0.987 

1.0017 

1.0024 
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theory  better  represents  experimental  shock  velocities,  with  Lagrangian  theory  yielding  too-low  values  of  V.  For  both 
clusters  of  experimental  shock  velocity  data,  3rd  order  Eulerian  theory  provides  a  better  fit  than  3rd  order  Lagrangian 
theory.  As  shown  in  Fig.  12(b),  all  4th  order  and  all  3rd  order  models  provide  a  reasonable  prediction  of  longitudinal 
shock  stress,  though  3rd  order  Lagrangian  theory  might  be  considered  overly  compliant  for  V/V0  <  0.95.  Second-order 
elastic  models  do  not  accurately  predict  Hugoniot  stress,  with  2nd  order  Eulerian  theory  too  stiff  and  2nd  order 
Lagrangian  theory  too  compliant. 

Considering  key  Findings  reported  in  Section  3.5  and  that  quartz,  sapphire,  and  diamond  all  have  B0>4,  the  lack  of 
definitively  greater  accuracy  of  Eulerian  theory  relative  to  Lagrangian  theory  is  somewhat  unexpected.  However,  quartz  and 
sapphire  are  not  cubic  and  have  directionally  dependent  covalent  and  ionic  bonding,  so  trends  from  the  ideal  model  of 
Section  3  that  enforces  Cauchy  relations  and  isotropic  symmetry  for  third-order  constants  may  not  apply.  Diamond  has 
strong  covalent  bonding  (unlike  the  noble  metals  for  which  anharmonic  properties  are  more  symmetric,  Hiki  and  Granato, 
1966),  and  examination  of  all  third-order  constants  in  Nielsen  (1986)  shows  that  (3.13)  and  (3.16)  are  not  well  respected. 
Nielsen  (1986)  reported  that  4th  order  Eulerian  theory  was  better  able  than  Lagrangian  theory  to  collectively  fit  atomic 
simulation  results  for  spherical  deformation  and  straining  along  [100],  [110],  and  [111]  in  diamond. 

Predictions  of  fourth-order  Lagrangian  and  Eulerian  theories  for  temperature  rise  0  (normalized  by  reference 

temperature  60),  entropy  jump  across  the  shock  r/  (normalized  by  specific  heat  at  constant  pressure  c§),  and  Hugoniot 

stress  P  (normalized  by  uniaxial  isentropic  stress  P)  are  listed  in  Table  5.  Isentropes  are  computed  as 

P1'  =  -J(Ci]E  +  J  CniE2  +  gCnn  E3),  (Lagrangian)  (4.60) 

Pn  =  -J-3(CiiD  + 1  CinD2  -(-  gCmiD3).  (Eulerian)  (4.61) 

Predicted  temperatures  are  similar  for  Lagrangian  and  Eulerian  theories,  with  temperature  rise  slightly  smaller  in  the  latter 
at  large  compression.  Predicted  entropy  production  is  positive  in  agreement  with  (4.9)  and  is  of  the  same  order  of 
magnitude  among  theories,  with  larger  y  predicted  by  Eulerian  theory  at  large  compression.  Recall  from  Section  4.1  that  the 
present  analytical  solutions  have  assumed  a  simple  form  of  specific  heat  wherein  the  contribution  to  internal  energy  from 
entropy  is  linear,  i.e.,  (4.31)  and  (4.43).  When  higher-order  Gruneisen  parameter  r\  =  0  and  c0kc^,  these  approximations  are 
most  accurate  for  t/<«2q).  From  Table  5,  such  conditions  hold  for  V/V0>0.92.  But  for  very  large  compression  (i.e., 
VyVo  =  0.88),  entropy  production  (especially  in  diamond)  is  large  enough  that  a  higher-order  representation  of  entropy, 
e.g.,  (4.21)  or  (4.22),  may  be  prudent.  Examination  of  stresses  in  Table  5  shows  that  P/P'  <  1.01  in  all  cases,  justifying 
isentropic  assumptions  used  in  previous  stress  analyses  (Fowles,  1967;  Clayton,  2009). 

Upon  examination  of  HEL  stresses  in  Table  4,  results  in  Table  5  are  deemed  valid  for  Z-cut  quartz  to  17/1^0  =  0.88. 
However,  for  X-  and  Y-cut  quartz  and  for  diamond,  the  HEL  is  exceeded  at  V /VQ  =  0.88,  and  results  are  most  valid  only  for 
17/170>0.92.  For  sapphire,  the  elastic  range  is  even  smaller  and  compression,  in  experiments,  is  elastic  only  for  17/170  =  0.96. 
Values  listed  in  Table  5  can  be  considered  extrapolations  when  compression  exceeds  the  HEL.  Above  the  HEL,  a  nonlinear 
theory  incorporating  dislocation  slip/twinning  (Clayton,  2009,  2010b;  Clayton  and  Knap,  2011a,  2011b)  and  cleavage  fracture 
(Clayton,  2006,  2010a,  2011b)  may  be  needed,  accounting  for  anisotropic  inelastic  deformation  mechanisms  and  their 
contributions  to  entropy  production  and  temperature.  If  large  densities  of  lattice  defects  are  generated  at  shock  pressures  at 
or  above  the  HEL,  consideration  of  their  effects  on  dilatation  (Clayton,  2009,  2011a;  Clayton  and  Bammann,  2009)  and 
tangent  elastic  moduli  (Clayton  and  Chung,  2006)  may  be  worthwhile. 


4.4.  An  alternative  lower-order  theory 

Purely  Lagrangian  and  Eulerian  theories  have  been  formulated  using  respective  strain  measures  E  and  D.  These  tensors 
are  not  the  only  possibilities  for  use  in  thermoelastic  potentials.  For  example,  a  generic  symmetric  strain,  in  material 


coordinates,  that  vanishes  when  F  =  1,  can  be  defined  as  a  linear  combination  of  E  and  D: 

G  =  (1~x)E  +  /D,  0<*<1.  (4.62) 

Lagrangian  (Eulerian)  theory  is  recovered  when  x  =  0  (x  =  1).  In  what  follows,  take  x=\<  giving5 

G='2(E  +  D)  =  \  (FtF-F“’F“t)  =  l(U2-U~2),  (4.63) 

where  U  =  RTF  is  the  right  stretch  tensor.  Free  and  internal  energy  densities  per  unit  reference  volume  are 

'F  =  W(G.  6),  U  =  U(G.rj).  (4.64) 

Thermodynamic  stress,  entropy,  and  temperature  obey 

S  =  d‘P /dG  =  dU /dG  =  S  ,  r}  =  —d'F/d6,  0  =  dU/di /.  (4.65) 


5  The  author  is  grateful  to  M.  Ortiz  for  suggesting,  before  calculations  in  the  present  work  were  undertaken,  the  possible  utility  of  a  strain  such  as 
(4.63). 
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Components  of  first  Piola-Kirchhoff  stress  are  obtained  from  chain  rule  differentiation: 


P  dU  dU  dGjj  U  r-  p-iF-iF-i 
P,tl  “  “  dC,j  dFkL  ~  2  S,;(5jiFw  +  FflI  F"t,Fl"i)’ 

where  0  can  also  be  replaced  with  'F.  Isentropic  second-order  elastic  moduli  are  defined,  in  Voigt  notation,  as 
Cl„  =  dl0/dGadGf. 


(4.66) 


(4.67) 


A  complete  three-dimensional  thermomechanical  theory  of  nonlinear  elasticity  can  be  formulated  in  terms  of  G 
paralleling  steps  taken  in  Section  2.  For  modeling  1-D  shock  compression,  a  scalar  theory  along  the  lines  of  those  in  terms  of 
E  and  D  in  Section  4.1  is  sufficient.  For  compression  parallel  to  the  X  =  X!  direction,  scalar  uniaxial  strain  is 

G  =  C„  =  1  (£„  +  D„ )  =  1  (F2-F~2)  =  IrV-D-  (4.68) 


For  the  present  purpose,  the  following  lower-order  thermodynamic  potentials  suffice: 
0  =  iC11G2  +  4o;?,  f  =  \CuG2. 


(4.69) 


In  this  representation,  temperature  is  constant,  thermoelastic  coupling  is  omitted,  and  isothermal  and  isentropic  second- 
order  elastic  constants  are  indistinguishable  and  are  both  represented  by  Cn.  Following  the  same  arguments  used  to  arrive 
at  (B.3),  this  second-order  elastic  constant  can  be  considered  equivalent  to  second-order  constants  of  Lagrangian  and 
Eulerian  theories  of  Section  4.1.  Hugoniot  stress  P  is  equivalent  to  the  second-order  elastic  isentrope,  and  is  calculated  using 
(4.66)  and  (4.68)  as 

p=-p11  =  -i(F+F-3)s=-ir3(i  +j4)du/dG=-ij-3n  +j4)gcu  =|r5d  +j4w-j4)cn.  (4.7o> 


This  relation  has  the  physically  appealing  feature  that  |P|->oo  as  J-»  0,  co. 

Predictions  of  (4.70)  are  labeled  “2nd  order  mixed”  in  Figs.  9,  11,  and  12.  Hugoniot  stresses  predicted  using  this  mixed 
Eulerian-Lagrangian  theory  are  more  accurate  than  those  predicted  using  purely  Eulerian  or  Lagrangian  3rd  order  models 
for  X-  and  Y-cut  quartz  [Fig.  9(a)  and  (b)[  and  X-  and  Z-cut  sapphire  [Fig.  11(a)  and  (b)[,  and  for  these  crystals  and 
orientations  are  of  comparable  accuracy  as  4th  order  theory.  Predictions  from  (4.70)  are  of  comparable  accuracy  as  those  of 
3rd  order  Eulerian  or  Lagrangian  approaches  for  Z-cut  quartz  and  diamond  [Figs.  9(c)  and  12],  but  in  these  two  cases  are  less 
accurate  than  4th  order  theory.  The  apparent  success  of  (4.70)  is  not  surprising  considering  that  2nd  order  Eulerian  and 
Lagrangian  models  tend  to  respectively  over-  and  under-predict  experimental  Hugoniot  stress  data.  As  noted  previously,  for 
many,  if  not  most,  kinds  of  single  crystals  only  second-order  elastic  constants  have  been  measured,  with  third-  and  higher- 
order  constants  unknown.  These  findings  suggest  that,  for  a  single  crystal  ceramic  or  mineral  undergoing  finite  elastic 
compression,  a  best  estimate  of  Hugoniot  stress  might  be  obtained  using  one-parameter  equation  (4.70)  if  only  second- 
order  constant  Cn  is  available. 


5.  Conclusions 

A  comprehensive  theory  of  Eulerian  thermoelasticity  has  been  formulated  and  compared  with  traditional  Lagrangian 
theory.  Analytical  solutions  have  been  compared  for  homogeneous  spherical  deformation,  uniaxial  strain,  and  simple  shear, 
for  cubic  crystals  with  fully  anisotropic  linear  properties  (three  independent  second-order  elastic  constants)  but 
directionally  independent  anharmonicity  (one  independent  third-order  constant).  For  a  typical  value  of  pressure  derivative 
of  the  bulk  modulus  of  four,  Eulerian  solutions  tend  to  demonstrate  more  physically  realistic  behavior  (i.e.,  rapidly 
increasing  stress  and  energy  at  very  large  deformation)  and  greater  stability  (e.g.,  positive  strain  energy  to  very  large 
deformation). 

An  analytical  solution,  accurate  to  fifth  order  in  strain,  has  been  derived  for  the  uniaxial  shock  response  in  Eulerian 
thermoelasticity,  paralleling  a  derivation  for  Lagrangian  theory.  Entropy  production  across  the  shock  is  minimally  third 
order  in  strain  for  either  theoiy.  Both  Eulerian  and  Lagrangian  treatments  provide  sufficient  accuracy  when  fourth-order 
elastic  constants  are  incorporated,  with  neither  theory  demonstrating  consistent  or  definitive  advantages  over  the  other  in 
describing  Hugoniot  data  for  quartz,  sapphire,  or  diamond.  A  second-order  model  incorporating  a  mixed  strain  tensor  that  is 
an  average  of  Lagrangian  and  Eulerian  strains  has  been  shown  to  provide  a  reasonable  approximation  of  the  Hugoniot  stress 
for  each  of  these  materials.  This  strain  measure  is  recommended  for  predicting  the  uniaxial  shock  response  of  other  hard 
anisotropic  crystalline  materials  if  higher-order  elastic  constants  are  unknown. 


Appendix  A.  Kinematics 

Let  x  and  X  denote  spatial  and  initial  material  coordinates  of  an  element  of  a  solid  body,  related  by  sufficiently  smooth, 
and  one-to-one  at  any  time  t>0,  functions 


x  =  x(X,  t),  X  =  X(x.  C). 


(A.l) 


J.D.  Clayton  /  J.  Mech.  Phys.  Solids  61  (2013)  1983-2014 


2011 


The  deformation  gradient  F(X,  t)  is 
F  =  V0X;  Fjj  =  dXi/dXj  =  djX  j. 

The  inverse  deformation  gradient  F-1(x,  t)  is 
F_1  =  VX;  Fj1  =  dX,/dXj  =  djX,. 


(A.2) 


(A.3) 


Partial  coordinate  differentiation  at  fixed  t  obeys  dj(  • )  =  F,ya,(  • ).  The  following  identities  (Clayton,  2011a)  that  follow  from 
FF-1  =  1  and  V0F  =  VqVqX  are  used  later: 


dFj/dFkL  =  -FTk'Fz/ 


djcFij  =  djFiK. 


(A.4) 


Current  volume  V  and  mass  density  p  are  related  to  initial  volume  V0  and  initial  density  p0  of  a  material  element  by  Jacobian 
determinant  /(X,  t): 

V 


=  —  =  J  =  detF. 
vo  P 

Lagrangian  Green  strain  £(X,  t)  is  defined  as 
E  =  \  (FtF-1  );  Eij  =  ^djXkdjXk-8jj). 

Letting  u(X,  t)  denote  displacement  and  Uj  =  UiSy, 

V0U  =  F— 1;  djllk  =  Fkj-8k]\  Ey  =  2(d!Uj  +  djU)  +  diUkdjllk). 
The  following  identities  apply  (Clayton,  2011a): 


§1  =  l<SL,FkJ  +  8IJFkl),  jL=JF-L{, 


17-1  17—1 
tJk  ■ 


dFkL  2  ii’  dpkL 
Eulerian  strain  D(x,  t)  is  defined  as 

D  =  j(l-F  ]F  T);  Djj  =  ^8lj-dkX]dkXj). 

Letting  u(x,  t )  denote  displacement, 

V U  —  I  F  ,  <)j  11  k  —  8jk  d,XK .  Djj  =  ^(djUj  +  djllj—  dkU[dkUj)8\i8jj. 

Comparing  (A.7)  and  (A.10),  and  noting  that 
djuk  =  diUk(8jj  +  djllj)  =  djUkSjj  +  0(||VU||2), 
it  follows  that  £  and  D  agree  to  first  order  in  displacement  gradients: 

D  =  E  +  0(||  Vu||). 

From  (A.4)  and  identity  <3  det  A/dAjj  =AJ)1  detA, 

dDij  1  , 
dFki 

From  definitions  (A.6)  and  (A.9),  for  ||£||  <  1  and  ||D||  <  1,  the  following  series  apply: 
D  =  1  [1-(1  +  2  £)-’]  =  l[l-(l-2£  +  4£2-8£3  +  •••)]  =  £-2£2  +  4£3— •, 

£  =  1  [(1-2D)-1  -1]  =  1[(1  +  2D  +  4 D2  +  8D3  +  1]  =D+2D2  +  4D3  + 

Consider  the  polar  decomposition  of  F: 

F=RU  =  VR,  RRt  =  1,  U  =  Ut,  V=Vt. 

Define  the  Eulerian  Almansi  strain  e(x,  t)  by 

e  =  1  (1  F  lF  1 );  e,j  =  \  (8ij-djXKdjXK)  =  i(a,  Uj  +  djUi-djUkdjUk). 


JJL  =  +  FJ'FJ3 ),  =JFk,FkJ. 


(A.5) 


(A.6) 


(A.7) 


(A.8) 


(A.9) 


(A.10) 


(All) 


(A.12) 


(A.13) 

(A.14) 
(A.15) 

(A.16) 

(A.17) 

When  there  is  no  rotation,  F  =  FT,  dtXj  =  dkXL8kjSiL,  and  dkUj  =  djuk,  leading  to  e  =  D  under  these  conditions.  From  the  polar 
decomposition  (A.16),  strain  tensors  can  be  expressed  in  terms  of  stretches  as 


£  =  l(l/2- 1),  D  =  l(l-ir2),  e  =  l(l-ir2). 

Spatial  velocity  gradient  l(x,  t)  and  its  symmetric  part  d(x,  t )  are 

l  =  Vu  =  FF“\  d  =  l(I  +  [T);  Ijj  =  djvj  =  djXj,  dy  =  f/djXj  +  djXj). 
It  follows  that 


(A.  18) 


(A.19) 


£  =  F  dF,  D  =  F1dF“ 


(A.20) 
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Of  particular  interest  are  situations  in  which  deformation  is  spherical  (isotropic): 

F=y1/3 1,  £“'  =J-1/31;  E  =  2  (/2/3-l  )1,  D  =  l(l-J“2/3)1.  (A.21 ) 

In  this  case,  the  following  limits  apply  as  J  =  V/V0-»0,  where  ||A||  =  (A  :  A)1/z: 

lim||£||  =  -v/3/2;  lim||D||->oo.  (A.22) 


Appendix  B.  Thermomechanical  derivations 


Material  coefficients  defined  as  derivatives  of  free  or  internal  energy  with  respect  to  £  are  related  to  those  with  respect  to 
D  as  follows.  First,  consider  thermal  stress  coefficients.  Setting  Fij  =  F'j =Sjj  in  the  reference  state,  (2.27)  yields 

Po*  =  ha,  foa  =  WJde)\o  =  Wjmo  =  P'o«-  (B.l) 

Recalling  that  c  =  c,  or  using  (2.36),  it  follows  that  Griineisen  parameters 

T0a=Ji0Jc0=f)0jc0  =  r0a,  T'0a  =  (dTa/dr,)\0  =  (dra/dn)\0  =  r'0a.  (B.2) 

Now  consider  second-order  elastic  stiffness  coefficients.  From  (2.52),  noting  also  that  P  =  0  in  the  reference  state, 

Co  a/i  —  Co  o/j>  C0„/;  =  C0a/).  (B.3) 

It  follows  that  for  higher-order  thermal  stress  coefficients  Poa/i*->ftoijKL’ 

PoijKL~PoijKL  =  _(dCjjK1/dd)|0  +  (dCyKL/d9)\0  =  P0IKSjL  +  PojKdtL  +  PoilSjK  +  Pojl^IK-  (B.4) 


Similarly,  for  higher-order  Griineisen  parameters  rQal!<^>r0ijKL, 

r  oijKL-r  oijkl  =  -&<)' (dClJKL/ dO)\o  +  9^  (dClJKL/ d6)\0  =  r0IKS]L  +  r  ojkSil  +  r0 ilSjk  +  r0 ilSjk-  (B.5) 

Note  that  by  definition,  T0a/!  =  [d(fla/c)/dEp]\0,  and  that  Ta^pap/c  in  general,  since  c  may  depend  on  £  when  p'a^0.  A  similar 
statement  applies  for  r0ali.  Finally  consider  third-order  elastic  coefficients.  Differentiating  (2.50)  with  respect  to  f  gives 
(Clayton,  2011a;  Clayton  and  Bammann,  2009) 


d3W 


nJkLmN  — 


cPPir  dAjjkL  d  _  _  —  — 

-  =  TF - (PiNpkM^NJML  +  ^ik^jO 

orki  dtmN  dtmN 


dFjj  dFkL  dFmN  dF  mN^F  kL  ormN  ormN 

=  F  ilF  kI<F  mM  C IJKLMN  +  ^imF  kK^JNKL  +  8 kmFilCjjNl  +  SikF  mM^JLMN 

=  FilF^F mM^ IJKLMN  “I"  & imF jq  Aq^ki  “1"  & kmF Nq^iJqL  “I"  &ikF jn  Am/Vni  —  (SimFjk  F NqPqL  +  ^km^Ni  FjJPqL  +  ^ikFjmFNgPqL)- 


(B.6) 


Differentiating  (2.51 )  with  respect  to  f  gives 

A-ijkLmN  =  ( Fpi  FRkFjr  ^Qr^Lt  F St  C PQRS~Fjk  PiL~FLj  P k] ~F jmF LmFkrtP In) 

or  mN 

F— 1  r—  1  u~  1  c—  1  £7—1  r-1  c— 1  £7—1  c— 1  A  sc—  1  £7—1  c— 1  £7—1  r— 1  c—  1  1  A  ,  £7—1  £7—1  c—  1  £7—1  c~  1  r—  1  E7—  1 A 

Pi  tRktUmtJr  tQrt^Lt  r st  P  NwtVw'^PQRSUV~(r  pm^  Ni  r  RktJr  FQrFlt  tS[  CpQ rs  +  rRmPNkPpi  f/r  ^St  c-PQXS 

' pF]t}lFNlrFplFRlFQ,rFljFsiCpQRs  +  Fq^F^F^F^Fj^F^F^Cpors 
im^Nt  Fpl  FRl  Fj ,r3  FQj  FSf  Cpqrs  4-  Fs^Fn J  Fp ’  FRl  F^  FqJ  F^1  C pqrs) 

~(Fjk  NlttiN  "F  F jj  A kJmN  "F  Fjq  F jq  FkpAipmN  +  SmkFjq  Fyq  Pjpj) 

+Fj^FNlPiL  +  Fp^FjqJPkj  +  fj/n Fpjq Fpg  FppPjp  +  FpmFnqFjqFkpP/p.  (B.7) 

Equating  (B.6)  and  (B.7)  produces  a  relationship  between  third-order  tangent  elastic  coefficients  at  any  deformation: 


Cr-— 1  £7 —  1  n—  1  c— 1  V~  1  c— 1  c— 1  u~  1  c— 1  n— 1  c— 1  ir—1  A 

IJKLMN  =  r,i  rpi  tjj  rQjtKktRktLl  tsl  r Tmt Nnt Un^ PQRSTU 

/  r— 1  r—  1  £7—1  £7—1  £7~1  £7—1  £7—1  £7—1  £7~1  £7—  1  A  1  £7~  1  £7—  1  £7~1  £7—  1  £7— 1  17—  1  £7—1  17—  1  £7—1  £7—  1  A 

Ii  t Ni  P Kk  P Rk  *  Mm** Pm ^ Jr  tQrtLt  tSt  ^ PQRS  +  P[i  Ppj  r Kk^ Mm** RmtJr  tQrtLt  t St  ^PQRS 
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This  equation  applies  analogously  for  either  isothermal  (C^kimn’^ukimn)  or  isentropic  (CyK1MN,CyKlMN)  coefficients.  In  the 
undeformed  stress  free  reference  state,  using  (2.50)  and  (B.3),  (B.8)  yields  the  following  relationship  between  third-order 
isothermal  elastic  constants: 


^ OIJKLMN  ~  Co  1JKLMN  ' 


SiKC 


+SJK^01LMN 

Similarly,  for  isentropic  constants, 


OJLMN 
djM  C  OINKL 


<5/lC 


OJKMN  ' 

■  SjlCoikmn 


■  SiMC, 


OKLJN 


+  SrnC 


■  SjNC0IMKL  - 


OKLJM 

SkmCi 


OIJLN  ' 


-  SKnC, 


01JLM  ' 


■  ^imC, 


0 IJKN 


-  sLNc 


0 1JKM- 


C  OIJKLMN  —  Cq  IIKLMN  +  SIkC 


'OIJKLMN  ' 


OJLMN  ' 


■  i5jiC" 


OJKMN  ' 


■  SiMC,' 


OKLJN  ' 


-  SjNC' 


OKLJM 


+sJkCqilmn  +  sJmC'oinkl  +  sJl^oikmn  +  sJn^oimkl  +  skm^oijln  +  sknCoIJLm  +  slmCqijkn  +  SLN^oijkm- 


(B.9) 


(B.10) 


These  derivations-which  effectively  equate  strain  energies  to  third  order  in  Taylor  series  with  respect  to  F-lead  to  relations 
among  third-order  elastic  coefficients  equivalent  to  those  derived  elsewhere  (Weaver,  1976;  Perrin  and  Delannoy,  1978)  by 
applying  series  approximations  [e.g.,  (A.14)]  and  equating  coefficients  of  like  terms  in  Lagrangian  and  Eulerian  energy 
potentials. 
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